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This thesis treats the subject of digital filters in two 
parts. First the analytical tools and synthesis techniques 
which are used for digital signal processors are explained 
and illustrated with examples. Emphasis is placed on 1) 
the generation and use of the digital transfer function, 2) 
design procedures in the frequency domain using methods of 
continuous-filter design as intermediate steps, and 3) 
real-time digital-filter implementation schemes. 

In the second part of the thesis a machine-language 
computer program is presented for the real-time simulation 
of digital filters on a hybrid computer. The program is 
described in detail and experimental results for several 
filter designs and realization schemes are reported. It is 
believed that such a computer program, with its inherent 
ability to accurately simulate a special-purpose computer and 
to provide external control of word length and clock frequen- 
cy, could be used to experimentally determine specifications 
for the LSI implementation of digital .filters . 
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I. 



INTRODUCTION 



The jyrocess of signal waveform filtering in linear sys- 
tems may generally be defined as the application of a linear 
operator -to an input waveform as it passes through a system. 
The transformation may occur in the time domain such as 
interpolation, extrapolation, differentiation and integra- 
tion; or -it may occur in the frequency domain such as low- 
pass filtering or smoothing, bandpass filtering, spectrum 
shaping and spectral decomposition. Operations such as 
these on continuous-time signals using analog filters are 
defined mathematically by linear constant-coefficient differ- 
ential equations. Analysis of these operations is made 
possible in the s domain by the LaPlace transform, and in the 
f domain .by - the Fourier transform. Similarly linear trans- 
formations^ on discrete-time signals may be defined by linear 
constant-coefficient difference equations, and discrete-time 
filters may-be analyzed and synthesized using the z-transform. 

A. A GENERAL DISCUSSION OF DIGITAL FILTERS 

A digital' filter is defined as a time-invariant linear 
operator on discrete-time signals. This linear operator is 
basically a computational process by which a sampled signal 
or sequence of numbers, called the input, is transformed into 
a second sequence of numbers, called the output. Occasionally 
the term "digital filter" includes the means for performing 
the analog-to-digital and digital-to-analog conversions. 
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However, since some signals are already digital in nature, 
"digital filter" in this paper will refer only to the compu- 
tational process described above. 

Digital ‘filters are time-invariant in the sense that they 
are described by linear constant-coefficient difference 
equations i Any given set of coefficients uniquely specifies 
a filter, and any change in these coefficients would neces- 
sarily result in a change in the specified filter. Thus, a 
digital filter is time-invariant. 

Implementation of the computational processes required in 
digital filters is accomplished using digital components such 
as shift registers, adders and gates. Several fundamental 
advantages to be gained using digital hardware are reduction 
in cost, higher speed, and increased reliability. Digital 
signal-processing techniques offer the following additional 
advantages over analog signal processing: 

1. Very predictable, drift-free performance of arbitrarily 
high precision can be obtained. This permits the 
realization of stable filters with very high Q's. 

2. There^are no impedance-matching problems since digital 
circuits are inherently isolated from each other. 

3. There is no restriction on critical filter frequencies. 
In particular, filters for frequencies of fractions of 
hertz are easy to construct. 

4. A greater variety of filter functions is possible since 
implementation is essentially a matter of numerical 
coefficients vice inductors and capacitors . 

5. Filter response can be changed easily by varying the 
proper coefficients. 

6. There is an intrinsic possibility of time-sharing the 
major filter sections. For example, a many-section 
filter^ can be implemented by building only one section 
and multiplexing that section by time-switching it. 

The filter coefficients can be changed for each section. 
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Alternatively, a given filter may operate on more than 
one input sequence by time multiplexing. 

These advantages in conjunction with the potential impact 
of large-scale circuit integration (LSI) on real-time digital 
filters assures wide application for these devices in such 
areas as audio systems, speech processing, seismology, ocean- 
ography, radar systems and data communications. 

B . HISTORICAL BACKGROUND 

The processing of discrete signals using linear filters 
or weighting sequences is traced back to the 1600 's by 
Kaiser [1] . He discusses the early use of classical niamerical 
methods arrd the subsequent evolution of the z-transform cal- 
culus. The theory of the z-transform and its application to 
the early digital signal processors are contained in the 
extensive literature on sampled-data control systems published 
in the 1950 's. See, for example. Refs. 1, 2, and 3. 

More recently the general-purpose digital computer became 
a vehicle for the simulation of continuous filters, and the 
development of digital-filtering techniques per se was greatly 
accelerated.' As computer technology advanced, particularly in 
the areas of memory size and computational speed, the real- 
time implementation of digital filters on general-purpose 
computers became possible.^ Now, with the advent of 



Indeed, the full impact of high-speed convolution via 
the fast Fourier transform (FFT) on real-time digital fil- 
tering is -not yet known. 
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large-scale circuit integration, the realization of digital 
filters in special-purpose hardware is becoming practical 
and economical. 

In short, the prospect of building real-time digital 
filters via LSI has generated considerable interest in the 
general field of digital signal processing. References 4-8 
contain much of the most recently published material on this 
relatively new subject. 
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II. MATHEMATICAL CONSIDERATIONS 



The analysis and design of digital filters requires a 
moderate degree of familiarization with the mathematics of 
discrete- time signals. In this section the theory of the 
sampling process and the z-transform calculus is reviewed, 
and their relation to the digital processing of signals is 
presented. Particular attention is paid to the methods of 
generating the digital transfer function, H(z); because, as 
will be shown later, this operation comprises the essence of 
the digital-filter design problem. 

A. THE SAMPLING PROCESS 

The operation of sampling a continuous-time signal can 
be thought of as a form of impulse modulation. If a signal 
f(t) is sampled every T seconds, then the sampling rate is 
'"Wg = 2'fr/T, and the sampled signal is a train of equally 

spaced impulses with varying amplitude given by 

00 

f*(t) = Z f (nT) 6 (t-nT) . (2-1) 

n=o 

In the analog-to-digital converter the amplitude of f(t) 
at each sample time, f (nT) , is "measured," and a digital 
word is then formed by encoding a sequence of binary digits 

t 

corresponding to that measured value. It is immediately seen 
that the acculracy of the "measurement" will be limited by the 
length of the computer word which the digital filter is 
designed bo handle. This phenomenon is known as quantization 
[4, 5], and its effects are discussed in Section IV, D. 
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Returning to Eq . (2-1), it is evident that, in the time 

domain, the sampled function f*(t) equals the original 
function multiplied by the modulating function, 6(t-nT) . 
Consequently, in the frequency domain, the spectrum of f*(t) 
will be equal to the convolution of an infinite number of 
equally spaced spectral lines with the original signal spec- 
trum. Hence, the periodic spectrum shown in Fig. 2-1 is 
obtained. It is important to note that in a real sampler the 
modulating pulses are not mathematical delta functions, but 
pulses with a finite width. This explains the decreasing 
amplitude in the spectra along the frequency axis on either 
side of zero frequency. The shape of this attenuating 
function is the familiar sin x/x shape, and its derivation 
can be formalized using the Fourier transform. 

Of particular significance is the phenomenon of fre- 
quency folding or aliasing. Again refer to Fig. 2-1. If 
w is the highest frequency component of the original signal 
spectrum and w < 2w , distortion will appear in the sampled 

o O 

signal spectrum because of the overlapping of signal compo- 
nents. In general, therefore, in order to recover the 
original spectrum with a low-pass filter, that spectrum must 
be bandlimited to ±to /2. This is known as Nyquist's Sampling 

o 

Theorem [1-3] and 2w is called the Nyquist frequency (minimum 
sampling frequency) . In practice there will always be some 
degree of frequency folding. Hence the total amount of in- 
formation is never recoverable. 
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Figure 2-1. Spectrum of a Uniformly-Sampled Signal 



It shoulti be emphasized at thisi^point that the sampling 
frequency must be the clock frequency for the digital filter. 
That is, all numerical computation that is required to gen- 
erate one output value must be completed in one sampling 
period. Hence, the computation time of the digital filter 
sets an upper limit on the sampling frequency. 



B. LINEAR DIFFERENCE EQUATIONS 

If an input sequence is de'noted by x , and an output 

' V ' . 

sequence by y^, then the digital filtering of to produce 
y^ may be expressed by linear constant-coefficient difference 
equations of the form 



Y + 
n 



bn y n 
1-^ n-1 






o n 



1 + 
1 n-1 



*"^^N^n-N* (2-2) 



Each term"in the above equation is a sequence of numbers. 

The sequence differs from y^ only in that it is delayed 

one unit (T) ^ in time. Equation (2-2) , written as a function 
of T, becomes 

(2-3) 

y(nT) + bj^y(nT-T) + •••+ bj^y(nT-MT) = a^x(nT) + a^x(nT-T) 

+ ••• + aj^x(nT-NT). 

Equations (2-2) and (2-3) are equivalent. The notation of 
Eq. (2-2) is often preferred because of simplicity. The 
interpretation of the generalized difference equation is as 
follows : ■ at time t=nT the new value in the output sequence 
y^ can be determined from the current input sample, the most 
recent value in the input sequence x^ , and a linear combina- 
tion of the past input and output values which are available 
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for computation. The constant coefficients, a^ and , are 

determined by the particular filtering function desired. The 

specification of these coefficients becomes the essence of 

the digital-filter design problem. 

An example of a first-order linear difference equation is 

\ 

y(nT) = b^y(nT-T) + x(nT). (2-4) 

(Here first-order refers to the fact that the highest order 

delay in the equation is one.) If x(t) is the complex 

sinusoid e^^^, then x(nT) is Accept, for the moment, 

“ i u)T 

that y(nT-T) can be expressed as y(nT)e ; i.e., y(nT-T) 
is y(nT) delayed one sampling period in time. (This idea 
will be formalized in the next section.) Then Eq. (2-4) can 
be written as 



y (nT) 



^ jno)T 




X (nT) e 

ej‘"T _ 




(2-5) 



Equation (2-5) shows that, in the steady state, the output 
is also a complex sinusoid, but modified by a "transfer 
function." The transfer function can be defined as H(e^^"^) . 
Its magnitude is 






(1 + b: 



- bj^coswT) 



1/2 



( 2 - 6 ) 



and the phase is 



'1' = COT - tan ^ ) 

cos(oT-b; 



(2-7) 



Equation (2-6) illustrates the frequency-selective properties 
of the transfer function, H(e^^'^) . Note that the magnitude 



. a 
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response is periodic in frequency with a period equal to 
OJg = 2ir/T. These results are, in fact, quite general and 
the concept of a discrete-signal transfer function will be 
discussed in detail in the next section. 

C. THE Z-TRANSFORM 

It has been stated that digital filtering involves the 
computational process of obtaining a solution to linear 
constant-coefficient difference equations of the form given 
in Eq. (2.2). Furthermore, it was implied that the steady- 
state solution to Eq. (2-2) for a sinusoidal input can be 
obtained by applying a transfer function to the input signal. 
Transformations which permit the formulation of a discrete- 
time signal transfer function are called sampled-data or 
z-transforms . The term "digital transfer function" will be 
used to describe a discrete-time transfer function of the 
complex variable z. 

1. Definition of the Z-Transform 

Given a discrete signal f, represented by a sequence 
of number-s f^^; its z-transform, F(z), is defined by the power 
series 

A " -n 

F(z) = Z , (2-8) 

n=o 

where z is interpreted as a complex transform variable. The 
variable z plays a role similar to that of the LaPlace 
transform variable s. The operation of taking the z-transform 
of a sequence is denoted by 

F(z) = Z(f^). (2-9) 
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As a simple example, if were the constant signal 

= 1, n = o, 1, 2, , (2-10) 

then the z-transform is found from Eq . (2-8) to be the 

geometric series 

00 _ 1 

Z(l) = Z = ZiJz|>|. (2-11) 

n=o 1 - z 

The theory of the z-transform is well documented 

[2,3]; detailed proofs and further examples will not be given 

here. Rathei?, an important theorem and several conclusions 

which apply to digital filters are the following: 

1. The z-transform of a delayed sequence can be derived 
from the definition. The result is 

: ‘I 

Z (f , ) = z"^Z (f ) . (2-12) 

n-k n ' 



2. The z-transform of the real part of a complex sinusoid 
is 



Z(a*^cos nb) 



1 V. -1 

1-a cos b z 



z|>a . (2-13) 



3. The z-transform of the imaginary part of a complex 
sinusoid is 



Z (a'^sin nb) = 



a sin b z 



-1 



(1-ae^^z ^) (1-ae 



(2-14) 



Note that the z-transforms of sinusoids yield poles 
in complex-conjugate pairs at a radius depending on 
the damping factor. 

4. The z-transform of any function of the general form 

f^ = n^a*^ sin(nb + c) (2-15) 



and linear combinations of these will always be a 
rational function of z“ , and such functions may be 
broken down into terms of the form of Eqs . (2-13) and 

(2-14) . This class of signals corresponds to the class 
of continuous-time signals with rational LaPlace 
trans forms . 
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The correspondence between the z-transform and the 
LaPlace transform will now be shown. From Eq . (2-1) a 

sampled signal is 

00 

f*(t) = ^ f (nT)6 (t-nT) . (2-16) 

n=o 



Taking the Laplace transform of both sides of Eq. (2-16) 
yie Ids 



F*(s) = 



n=o f(nT)e-^^" 






(2-17) 

This result is identical in form with the definition of the 
z-transform, Eq. (2-8), provided the substitution. 



-1 -Ts 

z = e , 



-1 



(2-18) 

is made. The complex variable z may be interpreted as 
a "unit delay operator" in Eq. (2-8) ; or, in a different 
sense , as an ordering variable whose exponent represents the 
position of a pulse in a sequence. This latter interpreta- 
tion, sometimes referred to as a- "generating function" [2] , 

h 

will become clear after the discussion of the z-transform 
inverse . 

2 . The Digital Tra,nsfer Function. 

The output of a digital filter is again written as 



Y + b,y , + •••+ b,.y = ax + a,x , 

-'n l-'n-l M-^n-M on 1 n-1 

+ ' ' * + a„x . 

N n-N 

Taking the z-transform of each term and applying the 
delayed sequence theorem gives 



(2-19) 
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( 2 - 20 ) 



-1 -M 

Z(Yn) + Z(y^) + ... + Z (y^) = a^Z (x^) 

+ . . . + 

Factoring the comition terms and simplifying Eq. 
result is 

? -i 
Z a • z 
1 



+ anz ^Z(x ) 
1 n 



a.,z ^Z(x ) 
N n' 



(2-20) , the 



YU) 

X(z) 



1 = 0 



1 + 



M 

t 

i=l 




H(z) - 



( 2 - 21 ) 



H(z) is defined as the ratio of the z-transformed output to 
the z-transformed input and will generally be called a 
digital transfer function. 

Another approach to the formulation of a digital 
transfer function is one using the notation of a linear 
operator H. If an output sequence is given by 



n 



= ax + a. X - + 
on 1 n-1 



( 2 - 22 ) 



which is simply the input sequence multiplied by a weighting 
function, then the notation can be simplified to 

y^ '■ = H (xj^‘) , k = n, n-1, • • * , (2-23) 

where H defines the filter. If h^ represents the impulse 
response of the filter, (i.e., the response to an input 
sequence which is unity at n=0 and is zero at all other 
times) , then the time response y^ is obtained from the 
discrete convolution theorem [4], 
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n % 

y(nT) = E x (kT) h (nT-kT) ; (2-24) 

k=o 

or, equivalently. 



y 



n 



E 

k=o 



X, h , 
k n-k 



(2-25) 



Taking the z-transform of both sides of Eq, (2-25) and again 
using the delayed-sequence theorem. 



Z(y ) = X Z(h ) + X, z ^Z (h ) + ••• , 

■'n o n 1 n ' 



or 

oo 

Y(z) = H(z) E Xj^z ^ = H(z)X(z) (2-26) 

k=o 

Again, H(z) is defined as the digital transfer function of 
the digital filter H. 

3 . The Z-Transform Inverse 

The transform of the output sequence Y(z) is seen to 
be obtained from the relation Y(z) =H(z) X(z) . Assuming that 
both H(z) and X(z) are known and are rational polynomials in 
z, there remains the problem of determining the sequence Y^ 
from Y(z) . Generally there are three methods for determining 
the inverse of a z-transform [2,3]. The first two are the 
use of a z-transform inversion formula and the partial frac- 
tion expansion. The inversion formula is the most general, 
the partial fraction expansion requiring factorization of 
the denominator polynomial into terms for which transform 
pairs are known. 

A third method, and usually the easiest, is the power 
series or long division method. Consider the function Y(z) 
written in terms of z ^ as 
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Y(z) = 



^ -1 ^ -2 ^ 
a + a^z + a«z + 
o 1 2 

+ b -z ^ + b' z ^ 
o 1 2 



(2-27) 



If the numerator is divided by the denominator using long 
division, the result is 



Y(z) = qQ + ^ + q^z ^ + *•• • (2-28) 

This power series gives the output directly. 

y*(t) = q^(t) + q^(t-T) + q2(t-2T) + ••• . (2-29) 

Here the variable z ^ is interpreted as the output sequence 
"generating function" mentioned earlier. Such a routine is 
easily performed on a computer or calculating machine. 

4 . Sinusoidal-Steady-State Solution of Linear Difference 
Equations 

In continuous-time linear systems the labor of taking 
inverse LaPlace or Fourier transforms can be avoided for the 
special case of sinusoidal input functions. The output of 
such a system in the steady state will always be a sinusoid 
of the same frequency, but whose magnitude and phase are 
modified by the transfer function of the system. The solu- 
tion for the magnitude of the output over a specified fre- 
quency range is called the frequency response of the system. 
This function is found by evaluating the system transfer 
function H(s) or F(jio) along the real frequency axis. It 
will now be shown that the response of a discrete-time linear 
system can be found by evaluating the digital transfer func- 
tion H(z) around the unit circle in the z-plane . Hence, 
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for the very common case of an input sequence being a 
uniformly sampled sine wave, the output sequence may be 
determined without calculating a z-transform inverse. 

In Section II, B Eq. (2-5) showed that the output of 
a discrete first-order system- to a sinusoidal input was the 

same sinusoid modified by a frequency-sensitive transfer 

. n (ji)T • sT 

function H(e-^ ) . Applying the substitution z=e , with 

s=jw, permits the transfer function to be written as a 

function of z. 

ioaT e^^^ 1 

> = tIt — - = ^ <2-30) 

e-^ - b, 1 - b, z 

1 1 



If the original difference equation, Eq . (2-4) , had been 

solved using the z-transform, the resultant digital transfer 

function would have been identical to the one in Eq. (2-30) 

above. These results are not just coincidental, but are 

required by the definition of the z-transform. Therefore, it 

can be stated generally that the frequency response of a 
2 

digital filter can be found by evaluating the digital trans- 

"1 (joT 

fer function H(z) of that filter at z=e-^ 

This discussion is illustrated graphically in Fig. 2-2. 
The mapping of the s-plane into the z-plane is accomplished 
by substituting s=o+jo) into Eq. (2-18). Then 



(a + io))T oT ju)T 

z = e =e •e-' , 



(2-31) 



2 

In this context the term digital filter is intended to 
refer only to frequency-selective filters in the analog 
sense, vice the general time-domain filters mentioned in 
Section I, A. 
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S-Plane (s=o+iaj) jto 





Figure 2-2. S-Plane to Z-Plane Transformation 
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where T=2 tt/uj For I jj I cu /2 the right porcion of the s- 
s s 

plane (positive o) maps into the area outside the unit circle 
in the z-plane. Similarly, the left portion of the s-plane 
(negative d) maps into the area inside the unit circle. Just 
as the frequency response of a continuous-time filter is 
determined by the value of its transfer function on the jco 
axis, the frequency response of a discrete-time filter is 
determined by the value of its transfer function on the unit 
circle. If the poles and zeros of a digital transfer func- 
tion are plotted in the z-plane, then the magnitude of the 
response at the frequency ujj in Fig. 2-2 is L 1 /L 2 L 3 and the 
phase IS ai-Uz-as. These results are obtained using conven- 
tional pole-zero graphical techniques. It should be noted 

that the magnitude response repeats after to reaches co , 2 w , 

s s 

etc., and is symmetrical around . This result agrees with 
the periodic spectrum of the sampled signal discussed previ- 
ously. In general the phase characteristic may or may not 
repeat [9]-. Also, the stability of a digital filter requires 
that the poles of the digital transfer function be located 
within the unit circle in the z-plane. 

The frequency response of a digital transfer function can 
be evaluated analytically by substituting e for z ^ in 

the functional expression H(z) . Assuming H(z) can be factored 
into products of second-order expressions, then 



^ Poles and zeros of a digital transfer function are 
found, as expected, by the factorization of the rational 
polynomials in z. 
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A + + A,e-^^“^ 

H(e-3“T, . 1 2 



1 _L D “j^^T . „ -j2o)T 

1 + Bj^e + 826 -’ 

A^e^“^ + Aj^ + Aje“2“^ 

e3“T ^ 



Expanding the exponential terms and regrouping yields 



_ - ..m [ (A,+A ) coscoT+A, ] +j [ (A -A-) sinwT] 

H(e ^ 2 i 2 ^ . (2-32) 

[ (B 2 +I) cosojT+Bj^) ] +j [ ( 1“B2) sinwT] 

Computer Program 1 may be used to evaluate the magnitude of 
the right side of Eq. (2-32) . The program is written in the 
FORTRAN IV language . 



D. S -PLANE TO Z -PLANE TRANSFORMATIONS 

Many of the digital-filter design techniques to be dis- 
cussed in Section III involve the digitalization of analog 
filters. These indirect design methods require transforma- 
tions from functions of s to functions of z. This section 
discusses the mechanics of various mapping transformations 
from the s-plane to the z-plane and compares the relative 
merits of each transformation. The formulas given for compu- 
tation of coefficients are not derived, but reproduced from 

the literature [6] . 

I k 

1. The Standard Z -Trans form 

This transformation is simply the extension of the 
definition of the z-transform. A proper rational function 
of s is expanded into partial fractions and the z-transform 
of each term is calculated. If RR represents the real part 
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of a residue in the partial fraction expansion, and RI 
represents the imaginary part, then real terms are trans- 
formed according to 



where 




1 + 



-1 



(2-33) 



A = T ( RR) 

o 

= - exp (uT) * 

Complex conjugate terms may be combined and transformed 
according to the transform pair 



RR t jRI RR - jRl ^ ^o ^1^ 

S - U - jv S - U jv 1 . t> “1 r. -2 ' 

1 + Bj^z + B 2 Z 

where 



(2-34) 



= 2T(RR) 

Aj^ = 2T exp(uT) [RR cos(vT)+RI sin(vT)] 

= 2 exp(uT) cos ( vT) 

B 2 = exp (2uTl . 

-1 -sT 

Since the transformation involves the relation, z = e , 
the mapping of the s-plane into the z-plane ta)ces place 
exactly as was previously described for Fig. 2-2. Notice 
that the transformation maps each strip in the s-plane, 
bounded by w == (n+l/ 2 )(j 0 g and o) = (n-l/2)a)^, where n is an 
integer, into the entire z-plane. Thus, the mapping of each 
succeeding strip is superimposed on the baseband mapping func- 
tion, If the transfer function H(s) has any significant 
response above o) = (i.e., if the product of the zero 

vectors divided by the product of the pole vectors is not 

negligible for w>a) /2, frequency folding will occur and the 

s 
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transformation may be unsatisfactory. Consequently the 
standard z-transform is normally used only for bandlimited 
functions. The problem is related to the one discussed for 
the sampled signal. In the present case, however, the 
transfer function itself must be bandlimited in addition to 
the input signal. In cases where H(z) is not sufficiently 
bandlimited a low-pass "guard" filter may be inserted in 
cascade to permit a successful transformation. 

2 . The Bilinear Z-Transform 

To circumvent the problem of frequency folding a 

transformation is sought which will map the entire s-plane 

into the z-plane. As a first step, the entire s-plane is 

mapped into an s' -plane bounded by the lines s' = ±jo3g/2. 

The mapping relation is 

s = I tanh (^) . (2-35) 

Note that this transformation also maps the entire s-plane 

into each succeeding strip bounded by co' = (n+H)o)'g and 

03 ' = (n-Js) 03 'g for integer values of n. This condition is 

necessary since the subsequent function of z must be periodic 

-s ' T 

in frequency. Then, substituting z = e into Eq. (2-35) , 

the bilinear z-transform becomes 



s 



2 

T 



tanh (^^) 



2 (1-e ®'^) _ 2 (1-z ^) 

— s * T — 1 * 

T (1+e ) T (1+z ) 



(2-36) 



This algebraic transformation does map the entire left half 
of the s-plane into the interior of the unit circle in the 
z-plane. Folding errors are eliminated since no folding 
occurs. The price paid for this feature is the nonlinear 
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warping of the frequency scale. 
Eq. (2-35) , then 



If s' = ju)' and s = juj in 



g)T 

2 



tan ( 




(2-37) 



Since co ' is the frequency which is really mapped into the 
z-plane, the warping of the frequency scale according to the 
tangent function results. However, satisfactory results can 
usually be obtained by prewarping the design frequency before 
performing the transformation. 

3 . The Matched Z -Trans form 

This transform generates a digital transfer function 
with poles and zeros matched to those of the continuous 
function. The mapping transformation is 

s e^'^ = z . (2-38) 

Then real poles or zeros transform according to 

s - Uj^ 1 + A^z ^ , (2-39) 

where 

= - exp(uT) 

Complex poles or zeros are combined into second-order factors 
and transformed according to 

(s-u)^ +’ >- 1 + A^z ^ ^ , (2-40) 

where 

- “2exp (liT) cos'( vT) 

A- = exp(2uT) 

The matched z- transform preserves the shape of the 
frequence response characteristics. Golden [6] points out 
that it may sometimes be necessary to multiply the resultant 
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“1 ri 

H(z) by terms of the form (1 + z ) to achieve satisfactory 
results. This operation is equivalent to inserting n addi- 
tional zeros at the half-sampling frequency into the digital 
transfer function. 

4 . Other Transformations 

Kaiser [1] discusses various other s-plane to z-plane 
transformations which may be used to achieve specific results. 
One particular transformation of interest is 

o 

2 - 2z COSO) T + 1 

s r ^ . (2-41) 

z"^ - 1 



This transformation, which is a variation of the bilinear 
z-transform, maps the entire imaginary axis of the s-plane 
onto both the top and bottom arcs of the unit circle in the 
z-plane. Setting s = 0 and solving the numerator of (2-41) 
for z shows that the origin of the s-plane is mapped into 
the points e ± Thus, applying this transformation to 

a low-pass H(s) yields a digital bandpass H(z) with center 
frequency of If we denote the analog cutoff frequency 

by and the digital cutoff frequency by then (2-41) 

can be written in equation form as 



j(i). 



j2w T jWpT 

e ° -2e costo T + 1 

o 



j 2 cjT , 

- 1 



or 






A 



COSO) T - cosa' T 
o D 



sinojj^T 



(2-42) 
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Again frequency warping is evident; but, as with the bilinear 
z-transform, prewarping may be used to produce useful results. 

5 . A Summary of S-Plane to Z-Plane Transformations 

Each of the transformations discussed in this section 
has advantages and disadvantages. In the design of digital 
filters the choice of a transformation must be dictated by 
these factors. 

First the standard z-transform, which is an expo- 
nential transform, requires the partial fraction expansion 
of the transfer function H(s) and the attendant factorization 
of the denominator polynomial. The transform does preserve 
the shape of the impulse time response; but, because of 
frequency folding, useful results can be obtained only for 
bandlimited functions of the low-pass and bandpass variety. 

The bilinear z-transform, being an algebraic trans- 
form, is easy to implement and requires no reduction of the 
transfer function in s. However, because of the frequency 
warping, this transformation should be used only for func- 
tions with piecewise-constant magnitude characteristics. The 
bilinear form does preserve flat gain characteristics. 

The matched z-transform, like the standard z- 
transform, is exponential in form, but requires only the 
factored form of H(s) , not the partial fraction expansion. 

It preserves the shape of the frequency response character- 
istics, and gives good results for high-pass and bandstop 
filters . 



/ 
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III. DIGITAL-FILTER DESIGN TECHNIQUES 



It has been shown that digital filters may be described 
by digital transfer functions of the form 



H(z) 



Y(z) 

XU) 



N -i 



E a^z 
i=o 



-1 



M 

1 + E b. z 
i=l ^ 



(3-1) 



Equation (3-1) is perfectly general and applies to all 
digital-filter functions. The nature of the filter function 
is defined completely by the a^^ and b^, and the specifica- 
tions of these coefficients becomes almost the entire filter- 
design problem. ' _ 

Digital filters are generally classified as being "recur- 
sive" or "non-recursive." Referring to Ea . (3-1) if at least 

one of the b^^ are non-zero, the filter is said to be of the 
recursive type; i.e., the computation of the present value of 
the output requires not only the present and past values of 
the input, but also the past values of the output (feedback) . 

If all the b. are zero, then the filter is said to be of the 
1 

non-recursive or transversal type; i.e., the output is simply 
a linear combination of weighted input values. This filter is 
also referred to as a "tapped delay line" filter. 

The design methods for the two classes of filters differ 
markedly and each class has its distinct properties. The 
non-recursive filter generally reauires a large number of 
terms to obtain relatively sharp transitions in frequency 
response. However, the non-recursive filters have excellent 
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phase characteristics, being quite linear ^in some cases. 

In contrast, the recursive filter can provide sharp cutoff 
characteristics with relatively few terms, but the phase 
characteristics are poor. 

In Section IV it will be shown that real-time implemen- 
tation of digital filters is presently much easier and more 
practical when the recursive filter structure is used. For 
this reason the discussion of non-recursive filters in this 
thesis is brief. All digital-filter design examples are 
recursive type filters to permit subsequent real-time 
implementation . 

A. NON -RECURS I \7E DIGITAL FILTERS 

The design methods for non-recursive filters fall into 
three categories. The first category is the classical 
numerical approach. In this method classical interpolation 
and differentiating formulas for equally spaced data can be 
digitalized directly by making correspondence between the 
unit advance operator "E" in numerical analysis and the 
z-transform operator. These filters are good only at low 
frequencies and they correspond to using the first few terms 
of a power-series expansion. This method is restricted pri- 
marily to filters whose function is integration, interpola- 
tion, differentiation, etc. 

The second design category involves the use of the Fourier 
series. If the magnitude of the desired frequency response 
is expanded in a Fourier series over |uj|<a)^/2 and the series 
is written in terms of cosines (or sines) the result is 
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H(w) 



OO 



(3-2) 



- ^ (nwT) . 

n=o 

„ . . -1 -jwT, Eg. (3-2) can be written as 

But, since z = e 

OO 

H(z) = a^ + j E 
"^n=l 

which is the required filter. The filter has an infinite 
number of terns and is realizable only if an approximation is 
made by truncating the series. The approximation is good 
for a rapidly converging series, but will be poor if there 
are discontinuities or steep transitions in the desired 
frequency response. 

The third category of design methods for non-recursive 
filters involve attempts to achieve a better convergent 
series, or to systematically weight error functions. Methods 
exist which utilize the Chebyshev series, least-squares 
approximation and "window" functions. Kaiser [1] treats this 
subject extensively and he presents several examples illus- 
trating the design methods. 

». RECURSIVE DIGITAL FILTERS 

Recursive digital-filter designs may be accomplished 
either in the time domain or in the frequency domain. His- 
torically designers of sampled-data filters were concerned 
with step-response, peak over-shoot, etc., and most transfer 
functions were analyzed and synthesized to satisfy time- 
domain criteria such as these. More recently, since the 
scope of digital signal processors has widened, more emphasis 
has been placed on steady-state response in the frequency 
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domain . Goid—and— Rader — ['4-}' have~forma-li' 2 e(J-v-an-ious— procte- 
du re s— f or- - d i g.i.r,a L- f i .1 tg r des iq n i n . th e_.fr eq ue n cy doma in- 
Golden- f 6 Whas— also._co ntrx buted. to-,the-developmeD±. these 

me-fehode-;- -'f-ur-thermorer' he'-has- written..-^. comp uter progr am 
whi-eh-'esserTt ±a’i*ly' pe r f or m&--a IJL- -the -de sign „,c,QjQj 3 .utati on s-» t o 
be--found"“in^ this“sectron’’. The^rx5g^am---wa&-Hiprtr~a:valrhab-ie---to 
the author^of" thl s^thesi-s*-: 

There are two general approaches to the design of recur- 
sive filters. They are the direct design method and the 
indirect design method. This latter method, which appears 
to be more popular, uses many of the standard continuous- 
time filters as intermediate designs. 

1 . Direct Design Method 

This method employs computational techniques to 
determine the coefficients, a^ and b^, directly from the 
filter specifications. Some methods involve the use of 
polynomial-approximation theory. Another technique is to 
locate poles and zeros directly in the z-plane , and then 
possibly to perform an iteration or trial-and-error proce- 
dure. These design methods might be employed where the 
specified response is not of a' standard form (low-pass, 

2 . Digital-Filter Design via Analog-Filter Design 
The indirect method of digital-filter design is 

accomplished in two steps. First an analog filter H(s) is 
designed to satisfy the filter specifications. Then a 
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transformation of the analog function is performed to d:eT:er- 
mine a digital transfer function H(z). In this manner use 
can be made of the many design and synthesis techniques for 
continuous filters which have been developed over the years. 
R o - fe r on OG ' ■ 1 0 — io ' ■ tho-^xt^n^-i ve-Hrre^rtmeTttr "of'^'thi s 

subgexxt^in, 

a. Topics in Continuous -Filter Design 

Standard- forras.^of-contrinuousfi'lters“have~-been’' 



developed— and- are^'weil 'documented ‘in- the literatufe 



The filters are referred to by the names Butterworth, 
Chebyshev, Elliptic, Bessel, Lerner and others. Each type 
of filter has certain characteristics which determine the 
most appropriate use for that filter. Of these many filter 
types the Butterworth and the Chebyshev designs are now 
briefly outlined to provide a basis for the digital-filter 
design examples which follow. 

CD The Butterworth Filter . The Butterworth 
filter can be specified by the relationship 



H(jco) 



1 + (w/w^) 



2n ' 






where to^ is the cutoff frequency. |H(jto)p is the squared 
magnitude of the filter transfer function. It can he shown 
that the transfer function H(s) is a rational function of s 
with a constant numerator and a denominator polynomial of 
degree n. For example, a fourth-order Butterworth filter 



with normalized to 1 rad/sec has a transfer function 



0 ) 




Figure 3-1. Response of a Butterworth Low-Pass Filter. 

Radian frequency , u) 




Figure 3-2. Response of a Chebyshev Low-Pass Filter, 
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H(S) = 



+ 2.613s^ + 3.414s^ + 2.613s + 1, 






The coefficients for other-order filters are readily found 
in tables Plots of the Butterworth frequency response 

for three values of n are shown in Fig. Generally the 

response is characterized by a flat passband and a stopband 
attenuation of 20 ndb/decade. 

(2) The Chebyshev Filter . The Chebyshev 
filter is specified by 

2 1 



H( ju)) 



1 + e^V^(io/io ) 
n ^ c 



Cl-i) 



where is a Chebyshev polynomial of order n. The Chebyshev 
filter is characterized by an equal ripple in the passband 
and a monotonic decay in the stopband. The ripple amplitude, 
6 , is given by 

1 



6 = 1 - 









The transfer function of a Chebyshev 
filter is also a rational function of s. The coefficients of 
the transfer function are determined by the specified value 
of e and are calculated by a somewhat involved but well 
established procedure A plot of the response of a 

typical third-order Chebyshev filter is shown in Fig. 
Generally the Chebyshev filter will have a steeper slope in 
the stopband than a Butterworth filter of the same order. 

(3) Frequency Scaling . Since most con- 
tinuous filters are designed using normalized frequencies 
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(i.e., = 1 rad/secj , frequency scaling required to 

shift the response of the filter up to the frequency band 
of interest. If we denote the desired cutoff frequency 
by then to shift the frequency response of the normalized 



transfer function to co make the substitution 

c 



n 






^-5 

(3-e*) 



(4) Frequency Transformations . High-pass, 
bandpass, and bandstop filters may all be obtained from a 
low-pass design by effecting the following frequency 
transformations : 



Low-pas s-to-high-pass 



to 



Low-pass -to-bandpass 
s ^ 



CO 

o ,s_ 
BW ^(0 

o 



lO 



C 



(^1 

rr=9b'^ 



Low-pas s-to-bands top 



BW 



to 

, s . o» 

CO ( — + ) 

O CO s 

o 




C3-9o4 





In the expressions ^ BW = oo „-co , and oo = / to « • to , 

C^CJ- o c ^ cx 

(the geometric mean) . 

These transformations are all performed on 

the normalized low-pass H(s) since they contain their own 

scaling. Note that when using the bandpass and bandstop 

transformations to is not necessarily in the center of the 

o 

band. Consequently final designs may be slightly asymmetric 
relative to the intended center frequency. 
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b. The Design Procedure 

A desired frequency response is normally spec- 
ified by the nature of the characteristic in the passband 
and the amount of attenuation in the stopband. For example, 
it may be required to design a bandpass filter which is flat 
in the passband centered at 1000 Hz, down 3 db at 800 Hz 
and 1200 Hz (cutoff frequencies) , and down at least 20 db at 
500 Hz and 1500 Hz. The flat characteristic would normally 
call for a Butterworth design, the 20-db attenuation spec- 
ification would determine the order of the filter and the 
3-db frequencies would specify the required frequency trans- 
formation. In the case where a final digital-filter design 
is required, the sampling frequency would be an additional 
and extremely important specification. 

The next decision to be made is the type of z- 
transformation to be used. Generally the standard z- 
transform would be used only for bandlimited functions to 
prevent frequency folding as discussed earlier. Bandlimited 
functions are typically functions of s whose denominator is 
of high degree and whose cutoff frequency is a small fraction 
of the half-sampling frequency. Usually the success of a 
standard z-transform design cannot be assured without 
actually testing the final filter. 

The bilinear z-transform works well when the 
frequency characteristic is relatively piecewise-constant . 
Hence, most Butterworth designs and small-ripple Chebyshev 
designs could be transformed satisfactorily. The algebraic 
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form of the bilinear 2 -transform makes it easy to implement; 
but prewarping of the critical frequencies in the s -domain 
is required to achieve the desired response in the z-domain. 

Finally the matched z-transform may be used to 
achieve a frequency characteristic that is high-pass, band- 
stop, or other than piecewise-constant . 

The digital-filter design procedures discussed 
above will be illustrated by examples in the next section. 

EXAMPLES OF RECURSIVE DIGITAL FILTERS 
The following examples were chosen to illustrate the 
digital-filter design methods which have been discussed. 

The order of the filters was kept low to facilitate the 
handling of the numbers. All calculations were performed 
on an electronic calculator, and the results were rounded 
to the sixth decimal place. Fe^f— mcTre~e^fflpJ.ex— fiJ.tejta_±he 
use_o£-~a-<;omput«r^~4Uio.g£OTVsuch:-sta:^tb^—J3ne™w£-2rfct-e-n-"by.- -Golden 
[ 6 ] ' 'WO u-ld"a-lmost-<3ejttaaJiiy^.a— •. 

EXAMPLE 1 . Design a low-pass digital filter to operate at a 
sampling frequency of 2000 Hz. The cut off (3-db) frequency 
is to be 200 Hz and the response is to be essentially flat in 
the passband. 

1. A second-order Butterworth filter is chosen with a sub- 
sequent transformation using the bilinear z-transform. 
Prewarping of the analog cutoff frequency is required. From 



44 






0) 



A 



= tan 



03 T 

(-^) 



L i- 



/ 



03^ = tan ( 2 tt • 200^- .5x10 ^ ^ (0.314159) = 0.324920. 

Note that the analog cutoff frequency 03^ has been prewarped 
to (0.324920) (2/T) = 207 Hz. 

2. A normalized second-order Butterworth filter is 

H(s) = . 

s + 1.414213s + 1 




( - 3-11 4 



Substituting^ s = s/0.324920 gives 

L 

0.105573 

H.(s) = — = . 

' s + 0.459506s + 0.105573, 

3. The bilinear z-transform is taken by substituting 

< 

s = (z-l)/(z+l) in Eq. -fflliiPy • The result is 






H^(z) = 



0.105573 (z^+2z+l) 



z^-2z+l+0. 459506 (z-1) (z+1) + 0.105573 (z^+2z+l) 



which reduces to 



H^(z) = 



0.067455 + 0.134910 z~^ + 0.067455 z ^ 
1 - 1.142980 z"^ + 0.412801 z”^ 






The frequency response of this digital filter was evaluated 
using Computer Program 1, and the results are plotted i«n 

The 3-db frequency is exactly 200 Hz. The poles 



The factor 2/T is dropped here and in the subsequent 
substitution for s to avoid handling large numbers. Should 
the factor be retained the result can always be reduced to 
the same final answer. 
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Figure 3-3. Response of a Low-Pass Digital Filter 



of the transfer function, C3BM) are: 

z = 0.571490± 30.293598 

These poles are well inside the unit circle indicating a 
stable filter. 

EXAMPLE 2 . Design a bandpass digital filter to operate at a 
sampling frequency of 2000 Hz. The cutoff frequencies, 
and 200 Hz and 500 Hz respectively, and the response 

is to be flat in the passband. 

1. A second-order Butterworth is again chosen to be used 
with the bandpass bilinear z-transform given in Eq. (2-41) . 
This transformation will yield a fourth-order function in 

z for a second-order function in s. 

2. Equation (2-42) is written again for convenience. 



6), 



cosw T - cosco_T 
o D 

sinWj^T 



1 ( 3 - 14 ) 



It is desirable to have the analog frequencies associated 
with each of the critical digital frequencies be negatives 
of each other. Then the digital transfer function which 
results by transforming H(s/|o)^|) will satisfy the response 
specifications at both critical frequencies, and • 

To accomplish this solve for cosoj^T in 



COSO) T - COSO)_,T 
o D1 

sinoJj^j^T 



- cos 



qoT 



+ cos 






sino)D2T 



(3-15) 



After considerable manipulation of Eq^ (3-15) using trigo- 
nometric identities, the result is 
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cos 



(3-16) 



COSO) T = 
o 



2 ^02'^ ^ 

cos i - 03j^2T) 



Equation (3-16) need be derived only once. It then becomes 
a formula for use in the design procedure for the bandpass 
bilinear z-transf orm, 

3. Returning to the example 



COSO) 



cos i (700 • ^ 1 T • T) 

T = -4 1 = 0 . 



509523 



cos (-300 • 2it 



T) 



(3-17) 



4. Then, using Eq. (3-14), 



-3, 



and 



t^DiT = (200) (2 tt) (0.5 -X 10 -^) = 0.628 32 , 



, 0.509523 - cos 0.62832 ^ cnocoo 

OJ = ' "r. /oooT ~ = 0.509523. 

A sin 0.62832 



(3-18) 



An identical result is obtained for ~ 1-57080 verifying 

the validity of Eq . (3-16) . 

5. Substituting s = s/0.509523 into the normalized second- 
order Butterworth transfer function gives 

H(s) = Q - • ■ ^59^14 . (3-19) 

s^ + 0.720574s + 0.259614^ ~ 

6. The digital transfer function is now obtained by substi- 
tuting s = (z^ - 1.019046Z + l)/(z^ - 1) into Eq . (3-19). 



The result after considerable algebra is 
H(z) = 



0 .131106-0 .262212Z ^-i-O . 131106 z“^ 



(3-20a) 



1-1.400064Z ^+1.272216z ^-0.658419z ^+0 , 272217z“'^ 



or 
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Figure 3-4. Response of a Bandpass Digital Filter 




(3-20b) 



(0,362085-0 .72 4170z~^-t-0 .362085z~^) ^ 

(l-1.222023z"^+0 ,603825z"^) (1-0 . 1780 41z"^+0 . 4 50 821z"^ ) 

The frequency response of the filter is evaluated using Eq. 

(3-20b) and the computer program. The results satisfy the 

specifications exactly and are plotted in Fig. 3-4. The 

poles of this bandpass digital filter are located at 

z = 0.611011 ± jO. 480094 
z = 0.089020 ± jO. 665505 

Again the filter is stable. 

EXAMPLE 3 » Design a bandpass digital filter with a frequency 
characteristic within the passband equal to the reciprocal of 
a gaussian or normal curve. (An application of this filter 
is given in Section VI, C.) Refer to Fig. 3-5.' The center ' 
frequency is 450 Hz and the sampling rate is 2000 Hz. 

1. The shape of Fig. 3.5b may be approximated by a second- 
order Chebyshev filter. (An nth-order Chebyshev filter has 
n peaks and n-1 troughs. Hence, a second-order filter with a 
large e will have a dip in the center of the passband.) The 
design of this particular Chebyshev filter will only be 
summarized here. 

A value of e was calculated to give the dip a shape which 
approximates the reciprocal of a normal curve. To do this a 
relation between the statistics of the curve and frequency 

5 

had to be established. The normalized filter was then 

^ In Section VI, C, the origin of the gaussian curve is 
seen to be the spectrum of a CW carrier which is frequency 
modulated by noise. Hence, the relation between the spread 
of the curve and frequency is determined by the actual 
hardware being used. 
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a. Gaussian Distribution in Frequency Domain 
y = mean = ce n ter -frequency 
a = standard deviation 




b. Reciprocal of Gaussian Distribution 




Figure 3-5. A Noise Filter Specification. 
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designed using standard methods [10]. Finally a frequency 
transformation was performed to give the required bandpass 
filter. The result is 

H(s) = 0.515422s^/(0.000000052s^ + 0 . 000039934s^ 

o (3-21) 

+ 1.250115878s + 283.7773716s + 2610409.654). 

The frequency response of this f ilter , ^ | H ( j w) j , is shown in 

Fig. 3-6. The results are quite satisfactory. The problem 

at hand is then to obtain a digital transfer function which 

will yield the same frequency characteristic. 

2. The matched z-transform is chosen to obtain the required 
digital transfer function. The poles of the transfer func- 
tion H (s) are 

s = -288.74512 ± j4639.3984 
s = - 95.23544 ± jl521.2546. 

Digital transfer function coefficients can then be evaluated 
using the matched z-transform relationships given in Eq. 
(2-40) . The results for the two second order transfer 
functions are 

ao = 1.0 

a^i^ = -1.0 

b^ = 1.178165 

= 0.749203 

and 

ao = 1.0 

a^^ = -1.0 

= -1381435 
= 0.909159. 



6 Figure 3-6 is the actual output from the Calahan net- 
work analysis computer program [11] . This program is an 
extremely useful aid in the present context of analog 
transfer function synthesis. 
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Response of a Digital Noise Filter. 



(3-22) 



Then the complete digital transfer function is 
H(z)= 



K(l-z d-z"^) 



(1+1.178165Z ^+0,749203z (1-1. 381435z ^+0.909159z ' 



where K is a gain factor. 

3. An examination of the frequency response of the transfer 
function given by Eq. (3-22) reveals a difference of 5.5 db 
in magnitude at the two peak frequencies. If a first-order 
zero at the half-sampling frequency is inserted in the 
transfer function, the desired symmetrical response can be 
better approximated. The response of this modified transfer 
function is plotted in Fig. 3-7 where the gain factor K was 
adjusted for 6— db attenuation at center frequency. 
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IV. IMPLEMENTATION OF DIGITAL FILTERS 



Let H(z) define an arbitrary digital transfer function 
derived by one of the methods discussed in Section III. The 
solution for the output transform, Y(z) = H(z)X(z), can be 
written as 



N _i ^ -1 

y(z) = E a^ z ^ X(z) - _ E b^ z Y(z) ; (4-1) 

i=o i=l 



or, in the time domain interpreting z ^ as the unit delay 
operator , 



n 



N 

= E 
i=o 



a . X 
1 n-i 



M 

1=1 



n-i 



(4-2) 



This linear difference equation can be realized directly 
using digital components. Refer to Fig. 4-1. Each rectangle 
inscribed with z ^ represents a one-sample delay. In real- 
time hardware this means a one-word memory whose capacity 
in bits is determined by the number of bits of accuracy used 
to represent and y^. Also a multiplier is required for 
each branch in the filter. In general these are variable- 
coefficient multipliers to permit changing of the transfer 
function. Notice that these multiplying coefficients are 
the same coefficients which appear in the transfer function 
H(z) . This is the basis for the earlier assertion that the 
design of digital filters is essentially complete once the 
coefficients have been specified. This is in sharp contrast 
to analog-filter design where i-he determination of the 
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Figure 4-1. A Recursive Digital Filter. 
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transfer function H(s) is often only a smal'l part of the 
design problem. 

The adders and multipliers together comprise the arith- 
metic unit of a special-purpose computer or digital filter. 
The implementation of Eq. (4-2) would involve clocking this 
computer at a rate f=l/T, and performing the indicated 
calculations during each sample interval. After 
the calculations are complete, the contents of the storage 
registers, or delay elements, must be shifted one location 
to be in the proper position for the next sampling period. 

For example, if the input value at time t = nT is stored 

as X , then this value must be stored in the location for 

n 

x^_j^ one period later at t = (nT+T) . 

In general, the expression for H(z) = Y(z)/X(z) is a 
ratio of polynomials in z ^ where the numerator is of degree 
N and the denominator is of degree M. If N > M then the 
computation of the current value of the output x^ in Eq. 

(4-2) would require knowledge of the input value For 

example, if N=5 and M=4, then the computation of y^ would re- 
quire the value of x^ among others. Since this is physically 
impossible, it can be stated generally that for a digital 
transfer function to be physically realizable the degree of 
the numerator must be less than or equal to the degree of 
the denominator. 

It is seen that the realization of digital filters in- 
volves basically a computational device to solve the linear 
difference equation which is defined by the digital transfer 
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function of the filter. However, the sequence of performing 
the arithmetic operations and the method of storing inter- 
mediate answers do have an effect on the efficiency of the 
filter and the accuracy of the results obtained.^ A digital 
filter is said to be in canonical form if, for a given order 
n, a minimiam number of adders, multipliers and delays are 
utilized. Three canonical filter structures are now 
discussed. 



A. THE DIRECT FORM 

Consider the general second-order digital transfer 
function 

, -1 -2 
a + a, z -t- a- z 

H(z) = ^—5 . (4-3) 

1 + bj^ z + b 2 z 

H(z) can be written as the product of two transfer functions, 
Hj^(z) and H 2 (z) . Then Eq. (4-3) becomes 



Y(z) _ 
X(z) 



1 

1 + b^ z ^ 4- b 2 






(4-4) 



and 

Y(z) = X(z)H^(^)H 2 (z) . 

Define the intermediate variable W(z) 
in the usual notation 



(4-5) 

X ( z) H^ (z) . Then, 



^ Efficiency refers to the speed of computation and the 
effective utilization of hardware. Both these items have 
become important "yardsticks" for measuring the usefulness 
of a digital filter. 
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Figure 4.2. A Recursive Digital Filter in Canonical Form. 
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(4-6) 



w 

n 



X 



n 



b 



1 



w 

n- 



1 



b„w 

2 n- 



2 






a w 
o n 



an w T 
L n-1 



+ a^w ^ 
2 n-2 



The general result is that the single difference equation of 
the form of Eq. (4-2) is converted mto two difference 
equations of the form 



w 



n 



X 



n 



M 

I 

1=1 



b . w 



1 n-i 



y 



n 



N 

I 

1=0 



a^w 



n-i 



(4-7) 



The schematic representation of Eq. (4-7) is shown in Fig. 

4-2. Notice that Fig. 4-1 and Fig. 4-2 are similar in form 
but Fig. 4-1 contains redundant delay elements. Hence, Fig. 
4-2 is a canonical filter structure. This means a more 
efficient use of hardware. 

Any nth-order digital filter may be implemented in the 
forms of Figs. 4-1 and 4-2. These are both called direct 
forms and are the simplest of the recursive filter structures. 
However, because of numerical considerations, the direct forms 
are normally satisfactory only for low-order filters. These 
numerical considerations will be discussed in more detail 
in Section IV, D. 



B. THE CASCADE FORM 

The general nth-order digital transfer function may be 
factored into a product of second-order transfer functions. 
The result is 
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H(z) 



(4-8) 



1 -1 -2 
m 1 1- a, 2 + a„ 2 

_ ^ li 2i 

” -1 -2 ' 
° 1=1 1 + z ^2i ^ 



where m is the integer part of (n+l)/2. Note that for H(z) 
with poles and zeros in complex-conjugate pairs the coeffi- 
cients and will always be real. The resulting 
difference equations are of the form 



Yj^(z) 


= H^(z)X(z) 


Y2(z) 


= H2(z)Yj^(z) 

• 

• 


Y(z) 


• 



by the same arguments which were made in the preceding section. 

The schematic representation of Eq. (4-9) is shown in Fig. 
4-^ It is seen to be a cascade of second-order subfilters 
and this structure has several important features. First 
the structure can easily be multiplexed. By using a basic 
second-order filter structure and appropriate timing and 
control circuitry/ a high-order filter can be implemented by 
time-sharing the second-order structure. The appropriate 
coefficients, and , are supplied for the m subfilters 
in time sequence. 

Also notice that if H(z) contains zeros on the unit 
circle, the coefficient in the subfilter which contains 

a second-order zero will be unity. This means that one less 
multiplication is required. Such a situation occurs fre- 
quently.' when the bilinear z-transform is used to design band- 
pass and bandstop filters. 
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Figure 4-3. The Cascade Form 




Figure 4-4. The Parallel Form 
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C. THE PARALLEL FORM 

The third canonical form of a digital filter is obtained 



by expanding H(z) 


in partial 


fractions of second-order 


terms . 


The result is 








H(z) = Yq + 


m Y , -r 


^li 


(4-10) 


1—1 1 3 T —1 , 3o —2 

liz + 2iz 


where y = a /b 
' o n' n 


Equation 


(4-10) yields difference 


equa- 


tions of the form 


Yo(z) = 


YoX(z) 






Yj^(z) 


(z)X(z) 






• 

• 




(4-11) 






H^(z)X(z) 






Y(z) 


m 

2 Y^(z); 

i=o 




and the filter is 


structured 


as shown in Fig. 4-4. The 



parallel form is comparable to the cascade form in perform- 
ance. The choice between the two might be made on the basis 
of the amount of algebra involved. Depending on the z- 
transformation employed the amount of work in arriving at a 
final filter structure can be considerable, especially for 
high-order filters. 



^Alternatively the expansion may be made of H(z) written 



as a function of z rather than z 



-1 
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D. NUMERICAL CONSIDERATIONS 



Until this point it has been assumed that both the cal- 
culation of digital filter coefficients and the manipulation 
of these coefficients in the filters themselves were done 
with infinite accuracy. Since this cannot be the case, and 
since even the input samples to the digital filter are only 
an approximation to the instantaneous values of the analog 
signal, the actual results may differ markedly from the 
theoretical results. 

As with any numerical process the accuracy of the cal- 
culated "answers" using digital filters is a very critical 
function of the hardware used and the computational technique 
employed. The finite computer word length is the basic 
limitation on accuracy. Errors introduced because of this 
limitation are: (1) finite precision in filter coefficients; 

(2) quantization of input data to a specific number of bits; 
and (3) round-off error in multiplications and additions. 
References 1 and 2 contain rigorous treatments of quantiza- 
tion effects in digital filters. The results are summarized 
below. 

Numerical problems can begin even with the digital filter 
design. The design methods of Section III involve factoriza- 
tion and partial fraction expansions. Generally it is 
recommended to perform these operations in the s domain and 
to digitalize the resulting subfilters. Although factoring 
of the polynomials in s is not required for the bilinear z- 
transform, better results are obtained if this is done. The 
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generation and manipulation of high-order polynomials in z 
should be avoided. In particular, finding roots of poly- 
nomials in z can be quite difficult because of the tight 
clustering of poles and zeros in the z-plane. 

The popular notion that a higher sampling frequency 
always produces better results is not true. Actually a higher 
sampling rate requires greater accuracy of computation. Per- 
formance may deteriorate if the sampling rate is increased 
beyond certain limits. 

Often the most critical numerical problem in digital 
filters is the quantization of the filter coefficients. 

These coefficients determine the locations of the poles and 
zeros for the filter and any perturbation in the coefficients 
causes a shift in the pole and zero locations. For high- 
order filters the poles are usually clustered in a small area 
near the unit circle in the z-plane. Perturbations in these 
coefficients may actually cause the filter to become unstable. 
Furthermore, the pole locations are much more sensitive to 
coefficient perturbations when using the direct filter 
structure. Consequently high-order filters are almost always 
realized in the cascade or parallel form. Weaver et al [7] 
present a procedure for determining the bounds on the sensi- 
tivity of pole positions to coefficient perturbation. Given 
these bounds the required coefficient accuracy can be cal- 
culated. 

The quantization of the input signal in the analog-to- 
digital converter is treated as a noise input to the filter. 
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The round-off error in the arithmetic unit is also treated 
as noise but its entry into the filter may occur in various 
places depending on the realization structure. In the direct 
and parallel forms the quantization noise due to round-off 
error is simply additive. However in the cascade form round- 
off error may cause unsatisfactory operation because of the 
multiplicative effect in a serial system. A complete anal- 
ysis of various implementation schemes can be accomplished 
using noise models such as those developed by Gold and Pader 
[ 4 ] . 

E. HARDWARE CONSIDERATIONS 

Any device which operates on signals in real-time accord- 
ing to a digital transfer function is, in essence, a special- 
purpose digital computer. The computer program, or 
equivalently the transfer function and implementation scheme , 
can be hardwired directly into the computer. The input data 
would then be the signal samples and the filter coefficients. 
It is anticipated that LSI will provide the mechanism for 
incorporating such a computer into a package that is small 
in size, economical in cost, and reliable in performance. 

The design of such a package is a formidable task, and several 
of the more important considerations are now discussed. 

The first consideration is the choice of the sampling or 
clock frequency. Unfortunately the requirements of a digital 
filter in this respect are self-opposing. To obtain sharp 
cutoffs, more narrow bandwidths , deeper notches, etc., 
higher-order filter functions are required with 
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correspondingly large computation times. Since the computa- 
tion must be accomplished in a single sampling interval, the 
sampling frequency must go down to provide a longer time T. 
However, as the sampling frequency is lowered, the usable 
frequency spectrum is reduced because the relation co <w /2 
must be maintained. 

The second major consideration is computer word length. 
The cost and complexity of a real-time digital filter varies 
in direct proportion to the number of bits which must be 
stored and manipulated. It appears that the key to success 
in reducing computer word length lies in developing new 
implementation schemes which are less sensitive to numerical 
inaccuracies . 

In general, therefore, it is of the utmost importance to 
design and organize special-purpose hardware which can per- 
form digital-filter type computations efficiently, accurate- 
ly, and as fast as possible. One of the more demanding 
computations in a digital filter is multiplication. In 
order to illustrate some of the problems in logical design 
peculiar to digital filters a binary serial multiplier is 
outlined in detail in Appendix A. The multiplier fbilows the 
general design proposed by Jackson et al in Ref. 8. 

F. THE FAST FOURIER TRANSFORM 

The methods discussed thus far in implementing digital 
filters all involve the technique of interpreting z~^ as the 
unit delay operator and obtaining a solution of linear 
difference equations directly in the time domain. This 
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technique is perfectly valid for analyzing the steady-st^te 
response of digital transfer functiohs to sinusoidal inputs, 
and has been shown to be very practical in real-time appli- 
cations. However, for complex filters, and in particular 
for non-recursive filters with a large number of terms , the 
number of multiplications and additions may’ be prohibitive. 

In these cases filter implementation by high-speed convo- 
lution using the fast Fourier transform (FFT) becomes 
desirable. 

Recall that the discrete convolution of an input sequence 
with a filter impulse response is given by 



y 



n 



M-1 



I h 
i=o 



n 



X 



n-i 



(4-12) 



where M is the number of samples being convolved. When M 
is large evaluation of Eq. (4-12) by z-transform methods 
becomes difficult and time consuming. 

Alternatively the sequence y^ may be evaluated by (1) 
taking the Discrete Fourier Transform (DFT) of the input 
sequence and the impulse response of the filter, (2) multi- 
plying these two transforms together, and (3) taking the 
inverse DFT of the result. This procedure is very similar 
to using the conventional Fourier transform to analyze 
continuous-time signals in the frequency domain. However 

9 

the theory of the DFT is somewhat more involved. 



This is due to the fact that the product of two DFT's 
the frequency domain is equivalent to circular convolu- 
tion in time, not linear convolution as intended in Eq. 
(^■^t2) . Nevertheless, the general results as stated above 
are correct. 
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The practicality of this method lies in the ability to » 

use the FFT techniques to evaluate the DFT ' s . The fast 

Fourier transforro is described in detail in Refs. 4 and 

12 . Briefly it makes use of the fact that under certain 

conditions the calculation of the coefficients can be done 

iteratively with a considerable saving in time. Generally 

the value M (the number of sample values) must be a highly 

composite number (contains many factors) . When M is a power 

of 2 the FFT requires a number of computations proportional 

2 

to Mlog 2 M vice M for the conventional convolution. 

Presently a large general-purpose computer is required 
to perform high-speed convolution via the FFT. It is antic- 
ipated that this algorithm will eventually be programmed into 
special-purpose hardware for real-time applications. 
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V. REAL-TIME SIMULATION OF DIGITAL FILTERS 



In recent years large general-purpose computers have been 
used extensively to perform the computations required in 
digital signal processors. Computer Program 1 is an example 
of the frequency-domain analysis of a digital transfer 
function. In the time domain the calculation of an output 
sequence of a digital filter can also be programmed for a 
general-purpose computer. In this case the actual difference 
equations are solved and the program becomes a simulation of 
the digital filter. One limitation in such a scheme is 
immediately evident; a large computer memory is required to 
store the input and output number sequences. Nevertheless, 
such computer simulations of digital filters do find wide 
application. Indeed, most algorithms involving the use of 
the fast Fourier transform are performed on a general-purpose 
computer . 

With the increased emphasis on the real-time implementa- 
tion of digital filters it becomes desirable to be able to 
simulate digital filters in real time. Such a simulation 
would provide an experimental basis for establishing specifi- 
cations for the realization of digital filters in special- 
purpose hardware. It would provide a means for investigating 
and verifying new design and implementation techniques. And 
finally it would provide a real signal to observe, measure, 
and analyze both quantitatively and qualitatively. 
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A real-time simulation of digital filters must contain 
certain features if it is to be used as an experimental anal- 
ysis and design tool. First it must be possible to control 
the two basic design parameters--clock frequency and computer 
word length. Precise control of the sequence of arithmetic 
operations is required to simulate different filter struc- 
tures. And finally the simulation scheme must contain the 
necessary input-output interfaces to permit the complete 
processing of a signal waveform in real time. In short, it 
is required that a simulation scheme be developed which can 
duplicate as closely as possible an arbitrarily designed 
special-purpose computer or digital filter. 

A real-time digital-filter simulation in the above context 
cannot be accomplished using only a large general-purpose 
computer. Aside from the obvious lack of any analog signal 
interface , the control of computer word length and precise 
arithmetic sequence is difficult, if not impossible, to 
achieve in most large-computer operational environments. 
Furtheimiore , the use of a large computer for relatively simple 
computational algorithms for an extended period of time is 
extremely inefficient and, hence, undesirable. 

On the other hand, the real-time simulation problem is well 
suited to a Scientifically oriented hybrid-computer system. 

Some of the characteristic features of such a system are: 

1. Built-in interface between analog and digital computers 
including analog-to-digital converters and digital-to- 
analog converters. 

2. Readily accessible analog functions such as amplifica- 
tion, scaling, clipping, integration, etc. 
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3 . 



Readily accessible digital functions such as clocking, 
counting, storing, calculating, programming, etc. 

4. Accessibility of such peripheral equipment as oscil- 
loscopes, vacuum-tube voltmeters, signal and function 
generators, line printers, card readers, tape recorders, 
and teletypewriters. 

5. Typical "hands on" operation of the computer system 
vice the "batch processing" type operation of a computer 
facility. 

In the following sections a real-time simulation for dig- 
ital filters is presented. The specific simulation scheme 
was designed for use on the hybrid- computer system located in 
the Department of Electrical Engineering at the Naval Post- 
graduate School. While the ensuing discussion must be 
specific in certain areas, the solution to the overall prob- 
lem is quite general. 



A. THE HYBRID-COMPUTER CONFIGURATION 

1 . The Physical System 

The hybrid-computer system at the Naval Postgraduate 
School consists of the following major subsystems: 

1. Xerox Data Systems XDS 9300 Computer 

2. COMCOR Ci-5000 Analog Computer 

3. Hybrid interface subsystem consisting of: 

a. Adage VMX multiplexer 

b. Adage VT-13 14-bit A/D converter 

Some of the characteristics of the XDS 9300 digital 
computer which are of specific interest are: 

1. 24-bit word including sign bit 

2. META-SYMBOL Symbolic Assembler 



This nomenclature was formerly 
terns SDS 9300 Computer." 



"Scientific Data Sys- 
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3. 1.75-microsecond cycle time 

4. 1.75-microsecond fixed-point add time 

5. 7.0-microsecond fixed-point multiply time 

6. direct memory access for A/D converter output 

7. integral 15-bit D/A converter 

8. 3 index registers available to the programmer 

9. priority interrupt system 

The Ci-5000 analog computer has a full compliment of 
operational amplifiers, potentiometers, passive analog com- 
ponents , logic components , clock sources , and computer 
control devices. These items will be explained as necessary 
in the discussion. 

The Adage A/D converter produces a 14-bit binary 
fraction including sign bit from an analog input signal with- 
in the range + 100 volts. This 14-bit fraction is stored in 
the most significant (leftmost) portion of an XDS 9300 24-bit 
memory location. That is, 

+100V = 37776000, 



8 



0 V 
-lOOV 



00000000 



40000000 



8 



8 



where the decimal point is assumed after the sign bit, and 
negative values are stored in two's complement form. 

2 . The General Simulation Scheme 

The digital-filter simulation which follows is built 
around a digital-computer program written in a symbolic 
programming language. This symbolic language, called META- 
SYMBOL, is the means whereby a machine-language program is 
generated in the XDS 9300 computer. 



} 
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A machine-language program is desirable for several 
reasons. First it provides the rigid computer control and 
accessibility that is required for real-time filter simu- 
lation. For example, arithmetic operations such as addition 
and multiplication are specified precisely in the desired 
sequence. The location in memory of each variable, machine 
instruction and data word is known, and the flow of digital 
information in the computer is completely specified by the 
program. Furthermore, a machine-language program permits 
the operator to monitor the execution, enter the program to 
insert modifications, or to stop the program to examine 
intermediate results. In general, the programmer not only 
tells the computer what to do, but how, where, and when to do 
it. This complete control of the situation is required for 
the accurate simulation of a special-purpose computer. 

The analog signal which is to be filtered is fed to 
the patchboard of the analog computer. It is amplified, if 
necessary, and then fed directly to a trunk line that is 
connected to one of the input channels of the A/D converter. 
Operation of the A/D converter is under control of the dig- 
ital-computer program. At a certain point in the program the 
analog signal is sampled and the quantized signal magnitude 
is stored directly in a previously specified location in 
memory. The digital computer then begins to perform the 
arithmetic operations of the digital filter. The filter, 
coefficients had previously been stored in memory and the 
sequency of operations is programmed to follow the filter 
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Figure 5-1. Block Diagram of the Hybrid-Computer Configuration 




implementation scheme exactly. Upon completion of the cal- 
culations the output value for that sample period is stored 
and also made available to the D/A converter. When the D/A 
conversion is complete the filtered analog signal is fed back 
to the analog patchboard where it is rescaled and applied to 
an oscilloscope or spectrum analyzer. Back in the digital 
computer the various input and output sample values which 
are used for computation are shifted in memory to be in the 
proper location for the next sample period. The computer 
then waits for a clock pulse to begin the next iteration of 
the entire sequence . 

The master clock is located in the analog computer. Each 
clock pulse causes an interrupt in the digital computer 
operation and immediately causes it to go to the starting 
point in the iterative program. Obviously the digital com- 
puter must have completed the calculations for a given sam- 
pling period before the next clock impulse occurs. Thus the 
computation time of the digital filter sets an upper limit 
on clock frequency. 

3 . The Hybrid-Computer Diagram 

A block diagram of the simulation scheme is given in 
Fig. 5-1. An audio signal generator is used as a signal 
source. Because of the +100V input range in the A/D con- 
verter the original signal is scaled by a factor of 10 in 
amplifier A034. This means the input signal must now be 
limited to + 10. Ov. The unsealed signal, also fed to 

the upper trace on the oscilloscope. 
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The output of the D/A converter, which is the filtered 

signal, is fed to a unity-gain amplifier A036 through poten-' 

tiometer P036. Assuming no additional scaling of e^^ was 

required in the digital computer a setting of 0.1 on P036 

should cancel the gain factor of 10 on the input side, and 

the magnitude of should be correct according to the 

transfer function of the filter. If a given value of e. 

produces overflow in the computer's arithmetic unit, the 

overflow indicator on the console is lit. In this case the 

signal input level can simply be reduced, or additional 

- 1-2 

scaling (usually 2 or 2 ) can easily be accomplished in 

the computer program. Potentiometer P036 is then adjusted 
to equate the net scale factor of the entire loop to unity. 

B. THE DIGITAL-FILTER SIMULATION PROGRAM 

Computer Program 2 is a general program which can be used 
for the real-time simulation of digital filters on the hy- ■ 
brid-computer system discussed in the preceding section. 

The source program is written in the META-SYMBOL language 
and appears on the right side of the page. The block of 
numbers on the left side of the page is the corresponding 
machine- language program generated by the META-SYMBOL assem- 
bler in the XDS 9300 computer. A complete description of 
the operation of this computer and the mechanics of the META- 
SYMBOL language will not be given here. Detailed information 
on both subjects is contained in Refs. 13 and 14 respectively. 
However, the following limited information is provided in 
order to explain the simulation program. 
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1 . Background Information on the XDS 9300 Computer 



The META-SYMBOL source program is originally punched 
on cards. Refer to Computer Program 2. Card column 1 is 
printed in the location of the asterisk in line number 1. 

The format for the symbolic program is as follows: 

1. Card columns 1-7 contain the label field. This field 
is used for address identification or labeling. 

2. Card column 8-15 contain the operation field. Either 
mnemonic machine instructions or assembly directions 
are entered in this field. 

3. Card columns 16-30 contain the operand field and may 
be used for various purposes. Put very simply, the 
operand answers the questions "what?" or "where?" 
regarding the operation field. 

4. Card columns 31-80 contain the comments field and are 
ignored by the computer. 

The computer word length in the XDS 9300 is 24 bits 
or eight octal digits. The rightmost eight digits in the 
block of numbers in the program printout represent the ma- 
chine instruction corresponding to the META-SYMBOL instruc- 
tion on the same line. Each machine instruction and data 
word in the program is stored in the memory location given 
by the first five digits in that line. For example, refer to 
line 34 in the program. The machine instruction 01600132 is 
stored in memory location 00035. Furthermore, the eight 
octal digits of the instruction word are formatted as follows: 

Digit 1 contains address and index register information. 

Digits 2 and 3 contain the instruction code. 

Digits 4 through 8 contain the address of the operand. 

An example will help clarify the discussion. Again referring 
to line number 34 in the program, the META-SYMBOL statement 
reads, "LDA ••• SUM+1." This means simply, "load the A 
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register with the contents of memory locatiqn SUM+1." The 
operation code for LDA is 16. Therefore, the machine in- 
struction states, "load the A register with the contents of 
memory location 00132." From line number 100 it is seen 
that the location in memory which is one greater than the 
location reserved for SUM is in fact location number 00132. 



addressing is somewhat more involved and will not be ex- 
plained here. Moreover, many of the instructions in the 
program are required simply to make the program work and do 
not warrant explanation beyond the comments in the printout 
itself. Most of the "subroutines" (ADOPS, etc.) are re- 
quired by the hybrid interface software. 

2 . Background Information on the META-SYMBOL Language 



involves the use of the computer ' s arithmetic unit to 
perform the prescribed calculations on certain sample values. 
The following list contains selected machine instructions 
from the META-SYMBOL program language which are used to 
simulate digital filters. 

Instruction Translation 



Address modification using indexing and indirect 



The most important part of the simulation program 



LDA 



Load the A register with the contents 
of the effective memory location. 



STA 



Store the contents of the A register in 
the effective memory location. 



ADM 



Add the contents of the A register to 
the effective memory location. 
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MUL 


Multiply the contents of the A register 
by the contents of the effective memory 
location. The full product is contained 
in the A and B registers. 


ALSD 


The contents of the A and B register, 
treated as a double-length register, are 
shifted left a specified number of bits. 


ARSA 


A right shift similar to 
involving the A register 


that above 
only . 


BRM 


Mark Place and branch. 




BRU 


Branch unconditionally. 




LDX 


Load the specified index 
a specified number. 


register with 


BRX 


Branch only if the associated index 
register has not reached its terminal 
state . 



The general use of the above instructions in the digital 
simulation are: 

1. LDA, STA, MUL, AND ADM are used to add or multiply and 
to store results in a desired location. 

2 . ALSD is equivalent to multiplying by 2^ where n is 
the number of bits shifted. 

3. ARSA is equivalent to dividing by 2*^. 

4. BRM is used to shift computer control to a desired 
"subroutine." The return to the program is automatic 
when execution of the subroutine is complete. 

5. BRU is equivalent to a FORTRAN GO TO statement and 
is used in this program to iterate a sequence of 
instructions indefinitely. 

6. LDX and BRX provide a "DO LOOP" mechanism whereby a 
certain group of instructions is performed a pre- 
scribed number of times. Indirect addressing is 
often used to shift an effective address for each 
new iteration. 

3. Fixed-Point Arithmetic 



The XDS 9300 computer performs either floating-point 
arithmetic or integer arithmetic. Since floating-point 
arithmetic is far more complicated, and since it would almost 
never be used in special-purpose hardware, integer arithmetic 
is the most practical alternative. However, in the digital 
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filter itself, fixed-point arithmetic is desired whereby 
a mentally placed decimal point in a computer word remains 
fixed throughout the calculations. Hence, a conversion 
scheme is required. 

.Recall that the basic word length in the XDS 9300 
computer is 24 bits including sign. The output of the A/D 
converter is a 14-bit word corresponding to a binary frac- 
tion of magnitude less than one. That is, the first bit is 
a sign bit and bits 2-14 are magnitude bits assuming a 
decimal point after the sign bit. The least significant 10 
■bits of a sample value in memory, therefore, are all zero. 
The word format maybe written symbolically as 

X . XXXXXXXXXXXXX 0 0 0 0 0 0 0 0 0 0 

where X indicates a one or zero. Since the magnitude of the 
input sample is less than one, this word can always be 
treated as a 23-bit integer (excluding the sign bit) in the 
computer, and the decimal point may be inserted mentally 
whenever it is desired to determine the real value. For 
example, if the fraction (0.5) = (0.01)2 is added to the 

fraction (0.25) = (0.01)2/ then the computer performs the 

addition 

010000000000000000000000 

001000000000000000000000 

011000000000000000000000 

20 

The answer is (11x2 - (6291456) . However if a decimal 

point is mentally reinserted after the sign bit, the correct 
answer is seen to be (0.110***0)2 = (0.75) If the 24-bit 

word in the answer above is input to the D/A converter, it 
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correctly interprets the value as 0.75 and the output is 75V. 
Thus, no conversion is necessary for data which is input via 
the A/D converter. 

The situation is quite different for the multiplying 

coefficients, a^ and . If a coefficient of (0 * 5 ) 2^0 were 

to be entered into the program as data, the integer (base 2 ) 

equivalent of which is zero, would be stored in 

memory. This is because the assembler converts all base 10 

data to base 2 integers when using integer arithmetic. The 

23 

solution is to multiply (0.5) by 2 . If this were done, 

then the 23-bit integer that is generated in the computer 
23 

will be 2 times too large which is exactly what is desired. 

(O.5x2^^)^0 = (4194304)^^0 = (010000000000000000000000)2 
Hence, if the base 10 integer 4194304 is entered as data the 
desired fraction is stored in memory. 

Since the a^^ and b^ for most cascade and parallel 
filter forms are within the range -2 to +2, a more convenient 
conversion is 

= c^x2^^ = c^(4194304), (5-1) 

where the general notation c^ includes all multiplying coef- 
ficients. The conversion in Eq . (5-1) is equivalent to di- 

viding the coefficient by two to preclude any overflow in 
the A register. After the multiplication is complete the 
answer is multiplied by two to obtain the correct numerical 
result. In Eq . (5-1) is the base 10 integer entered in 

the computer program. 



83 



4 . A Discussion of the Program 

Computer Program 2 contains the scaled coefficients 
for the second-order low-pass filter which was designed in 
Section Illf C. The digital transfer function for that 
filter is 



H(z) = 



0.067455 + 0.134910Z ^ + 0.067455z ^ 
1 - 1.142980z“^ + 0.412801z“^ 



(5-2) 



and the corresponding difference equation is 

®2='n-2 ‘ ‘'ly.i-l ‘ '"2yn-2- < = 



If the coefficients are scaled according to Eq. (5-1) 
the result is 



a^ = 0.067455 
a^ = 0.134910 
a2 = 0.067455 
=-1.142980 
b2 = 0.412801 



Aq = 282927 
= 565854 
A 2 = 282927 
= 4794006 
B 2 =-1731413, 



where the signs of the are reversed to permit addition 
instead of subtraction.^^ 

The sample values are stored in memory in the fol- 
lowing order 



DATA 


= X 

n 


DATA+1 


= X T 

n-1 


DATA+2 


= 0 

n-2 


SUM 


>1 

I) 


SUM+1 


1 — 1 

1 

c 

>1 

II 


SUM+2 


= ^n-2* 


of the 


signs c 



1 , 



applies to all implementation schemes and will be done for 
all remaining examples without comment. 
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Then the sequence of instructions from line 22 to line 41 
in the program solves the following equation: 



an Al A2 

SUM = • (DATA) *2 + (DATA+1) *2 + (DATA+2) *2 



(5-4) 

+ |i- (SUM+1) *2 + 1^* (SUM+2) -2 , 

which is equivalent to Eq. (5-3). Thus, Computer Program 2 
is a simulation of the first of the two direct filter struc- 
tures discussed in Section IV, A. 

Notice that all the coefficients were divided by two 
to permit a more general program. Also recall that, because 
integer arithmetic is being used, the most significant portion 
of a product is contained in the A register. (The entire 
product is in the A and B 48-bit register.) If the contents 
of the B register are dropped completely, a correct 24-bit 
product is still retained since the value is interpreted as 
a binary fraction. 

At the end of the arithmetic sequence the variables 
are shifted in memory to prepare for the next sampling inter- 
val. This is accomplished, by the instructions in lines 4 3 
through 46. This sequence is performed five times (until the 
contents of the index register become negative) , and the 
effective addresses of the load and store instructions are 
decremented by one for each iteration. The resultant shift 
in memory is shown below: 
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Shift Order of Occurrence 



DATA 


last 


y 




DATA+l 


5 th 


1 

DATA+2 


4 th 


y 




SUM 


3rd 


4- 

SUM+1 


2nd 


4 

SUM+2 


1st 



The current DATA+2 and SUM+2 values are lost which is of no 
consequence. The contents of SUM become correct in the next 
sample period after the STA instruction is used in line 25. 

Implementation of more complicated filters will be 
illustrated in the examples which follow. The only differ- 
ences in the program are the sequence of arithmetic opera- 
tions and the arrangement of the data and the coefficients 
in memory. 

In the beginning of this section on digital-filter 
simulation it was mentioned that it would be desirable to 
control the computer word length. This can be done conven- 
iently in the program by the use of the ETR (extract) 
instruction. Refer to line 64 of Computer Program 2. The 
ETR instruction performs a logical AND operation between the 
contents of the A register and the operand of the ETR in- 
struction. In the D/A subroutine the output value y^ is 
loaded into the A register. Then a logical AND is performed 
between y^ and the contents of MASK which are 

( 111111111111111000000000) 2 
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The result is that y is truncated to 15 bits which is 
the word size used in the D/A converter. If an ETR instruc- 
tion is inserted after each multiplication, then all the 
arithmetic in the digital filter is truncated to the number 
of bits specified in location MASK. This means that the 
effective word length of the computer can be controlled by 
altering the contents of location MASK. Additionally the 
accuracy of the multiplying coefficients can be reduced by 
truncating the number of significant bits in the storage 
locations. This would normally be done manually since the 
coefficients are constant in time. 
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VI . EXPERIMENTAL RESULTS OF REAL-TIME DIGITAL-FILTER 
SIMULATIONS 



The digital-filter simulation scheme described in the 
previous section was used successfully for the real-time 
simulation of various digital filters. The same basic pro- 
gram may be used for all filters. Different filter struc- 
tures are simulated by changing the sequence of arithmetic 
operations. Different filters of the same structure may be 
simulated by changing the filter coefficients in the program. 
For the remaining examples only the arithmetic and data- 
shifting sequences and the filter coefficients are shown in 
the referenced computer programs. 

A. A LOW-PASS FILTER 

A second-order low-pass filter with a cutoff frequency 
equal to 200 Hz and a sampling frequency equal to 2000 Hz was 
designed in Section III, C. The transfer function of that 
filter is 

0.067455 + 0.134910z"^ + 0.067455z"^ 

H(z) = = • (6-1) 

1 - 1.142980z“-^ + 0.412801Z 

This is the same filter that was used to describe the general 
simulation program. Computer Program 3 shows the arithmetic 
and data-shif ting instructions necessary to implement this 
filter in the second-order canonical form- Notice that in 
this case additional scaling of the input data became 
necessary. This was easily accomplished using the AREA in- 
struction in line 23. The difference equations which are 



# 




f = 200 Hz 
s 

Time = 2 msec/div 
Amplitude = 5V/div 



Figure 6-1. An Output Signal from a Low-Pass Digital Filter. 
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implemented by Computer Program 3 are 



SUM _ SUM , B1 ,SUM+1, ^ B2 ,SUMt2, ^ 

4 " 4 2*^ 4 ^*^■^2*^ 4 ^ 



Y AO ,SUM, ^ A1 ,SUM-rl. ^ A2 ,SUMi-2 
4 = — ^*2 - — *^-4 ^*2 ^ 4— 



(6-2b) 



and the data is stored a.n memory as 



SUM 


- w 

n 


SUM+1 


w T 
n-1 


SUM+2 


^n-2 


Y 


= ^n ' 



where w is an intermediate computational variable. 

Equations (6-2) are identical in form to the general 
second-order difference equations for the canonical form, 

Eqs. (4-7) . The first term on the right side of Eq. (6-2a) 
is labeled SUM only for convenience in the program. This 
location contains the current input value at the time line 
28 is executed. Also since Y/4 is the answer obtained from 
the program, P036 in Fig. 5-1 must be set to 0.4 to complete 
the magnitude scaling of the output. 

The actual input and output waveforms for this filter are 
shown in Fig. 6- la. The ripple in the output wave is caused 
by the presence of the sampling frequency and its harmonics. 
These high frequencies could be eliminated by additional 
analog filtering if a smooth waveform is required. Figure 
6-lb shows the attenuation of the output at 200 Hz. The 
complete frequency response of the filter is exactly identical 
to the frequency response calculated for the transfer function 
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fg = 200 Hz 




2000 Hz 



1000 Hz 



Figure 6-2. The Spectrum of the Output Signal. 
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That response was plotted ±n Fxg . 3-3. The phase shift, 
caused by the insertion of tne filter is also evident in 
Fig. 6-1. 

Figure 6-2 shows the spectrum of the output waveform for 
an input signal frequency of 200 Hz ; The presence of the 
spectral component at 1800 Hz (2000 Hz - 200 Hz) verifies 
that the frequency spectrum of the filtered signal is indeed 
periodic with a period equal to the sampling frequency. 



B. A BANDPASS FILTER 

A bandpass digital filter was designed in Section III, C. 
The filter characteristics included a flat passband from 
200 Hz to 500 Hz for a sampling frequency of 2000 Hz. The 
transfer function of the filter is 

,,, , 0. 131106-0. 262212z“^+0.131106z"^ 

H(z) = ^ 

1-1.400064Z +1.272216Z -0.658419z -^t0.272217z 



1 . Simulation of the Direct Form 

Computer Program 4 gives the abbreviated coding for 
a direct form of this filter. The difference equation, in- 
cluding the in-program scaling of line 23 , 



SUM 

4 



AO ,DATA, A1 ,DATAt1, ^ 

2M4 ) 'Z + 2*^ 4 '*^ 



A2 .DATA 1-2. 

—•‘—4 ' '2 



^ M. (DATAil) .2 . M. ,DAT|^, . 2 . -2 



(6-4) 



Figure 6-3 shows the output waveforms of this filter at 
several frequencies. The complete frequency characteristic 
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fg = 200 Hz 
Time = 1 msec/div 
Amplitude = 5V/div 



fg = 320 Hz 
Time = 1 msec/div 
Amplitude = 5V/div 



Figure 6-3. Several Output Signals from a Bandpass 
Digital Filter. 



fg = 650 Hz 
Time = 1 msec/div 
Amplitude = 5V/div 
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does agree with the plot of the transfer-function character- 
istic given in Fig^ 3-4^ 

2 . Simulation of the Cascade Form 

The cascade filter form requires the factorization of 
the digital transfer function into a product of second-order 
subfilters. The result for the present example is 



H { z) = 



( 0.3 620 85-0, 724 170 z"^+ 0.36 20 85z ^ 



(6-5) 



(l-1.222023z~^+0,603825z ^) U-0 , 178041z"^+0 .450821z 



The implementation of the cascade form involves the 
routing of an input signal through a cascade of subfilters 
(two in this case). If the same basic second-order filter 
is used for all the subfilters / then when the signal comes 
out of one subfilter it must be routed around to the input 
of the common filter as the filter coefficients are changed. 
This is the most appealing feature of the cascade form. 
High-order filters can be implemented by time multiplexing a 
low-order filter section. 

The real-time simulation scheme can be programmed to do. 
the same thing. Refer to Computer Program 5. A basic 
second-order canonical- form filter is coded with an addition- 
al indirect addressing feature. Index register 3 adds either 
one or zero (depending on which of the two subfilter iter- 
ations is in progress) to all operands which are followed by 
the characters ",3", The data is stored in memory as follows: 

SUM = output of first subfilter = input to second subfilter 
SUM = input sample = input to first subfilter 
SUM+1 = delayed sum of second subfilter 
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SUM+1 = delayed svam of first subfilter 
SUM+2 = delayed sum of second subfilter 
SUM+2 = delayed sum of first subfilter 

Y = final output of complete filter 

Y = output of first subfilter, 

and the coefficients are stored as shown in the program. The 
sequence of instructions from "REPEAT" to "BRX REPEAT" are 
the basic arithmetic instructions for a second-order canon- 
ical filter structure. This sequence is executed twice 
during each sample period. Indirect addressing causes the 
proper storage locations and filter coefficients to be used 
for each subfilter. In this way a time -multiplexed filter is 
simulated. Care must be taken to insure the proper shifting 
of the data between sample intervals. Refer to lines 51 

through 56 in the program. 

The bandpass digital filter under consideration was 
simulated using Computer Program 5. The resultant frequency 
response was identical to both the calculated response and 
the response of the direct-form simulation. 

3 . Simulation of the Parallel Form 

The parallel form of a digital filter may be chosen 
if the paramount consideration is speed in lieu of conserva- 
tion of hardware. In this scheme the total response is the 
sum of the responses to the individual subfilters which are 
arranged in parallel. Refer to Fig. 4-4. Evidently the cal- 
culations for each subfilter can be performed simultaneously 
if no multiplexing of hardware is desired. 
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The implementation of the parallel-niter structure 



requires the partial fraction expansion of the digital trans- 
fer function into a summation of second-order transfer 
functions. This operation may involve considerable labor if 
the order of the filter is high. Computer Program 6 is a 
relatively simple FORTRAN program for obtaining the partial 
fraction expansion of some rational polynomials. The multi- 
plicity of poles of the rational function cannot exceed two. 
Fortunately digital transfer functions normally do not have 
multiple poles, so the program can be used. Minor modifica- 
tions to the program would permit expansions of arbitrarily 
high-order transfer functions. 

The bandpass filter presently being considered is of 
fourth order. Accordingly, the transfer function given in 
Eq. (6-3) must be expanded into the form 



a, ( z) a, ( z) 
H(z) = y + ^ ; ■ + 



6^ (z) 



^ 2 ) 



( 6 - 6 ) 



The algebra will not be given here. The results are: 

Y ^ 0.131106 

a^_(z) = -0.383754Z -t- 0.149628 

B^(z) = z^-1.222023z + 0.603825 

OL^(z) = 0.567310Z -I- 0.046308 

82 ( 2 ) = z^-0,178041z -t- 0,450821. 

Notice that the expansion was made for H{z) equal to a func- 
tion of z rather than z Then, after Eq. (6-6) is written 

in standard form (a function of z ^) , the filter and program 
coefficients become 
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Y 



0.131106 



GAMMA = 549898 




-0.383754 



0.567310 

0.149628 



Al^ = -1609581 
2 = 2379471 



0.046308 

-1.222023 



A2j_ = 627585 

2 = 194230 



-0.178041 

0.603825 

0.450821 



Bl^ = 5125536 

2 = 746758 



B2^ = -2532626 



= -1890880 



Computer Program 7 is the simulation of the parallel form 
of the bandpass digital filter. The program is identical to 
the program for the cascade form except that the outputs of 
each individual subfilter are stored separately and summed 
together at the end. The simultaneous operation of the 
subfilters obviously cannot be simulated because the com- 
puter's single arithmetic unit must be time-multiplexed . The 
correct results for the parallel filter are obtained although 
the increase in operating speed is not. In fact, the parallel 
filter simulation takes slightly longer than the cascade 
simulation because of this additional summation sequence 
(lines 49 through 53) which is required at the end. 

The results of the real-time simulation of the bandpass 
filter using Computer Program 7 were essentially identical 
to those of the previous two simulations. 

4 . Experimental Observations of Quantization Effects 



decimal-place accuracy in the design of the bandpass digital 
filter. This accuracy is retained by the 23-bit magnitude 
representation in the XDS 9300 computer. If the least 



The multiplying coefficients were calculated to six- 
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significant bits of the coefficient data words are system- 
atically truncated, some insight into the effects of 
coefficient quantization can be obtained, A software package 
is available with the XDS 9300 computer which permits alter- 
ing the contents of memory locations via manual input to the 
teletypewriter. 

The coefficients for each of these bandpass-filter 
configurations were truncated in steps of single octal 
integers with the proper rounding being performed. Generally 
the effect of the less accurate coefficients was to gradually 
alter the frequency-response characteristic unril it no 
longer satisfied the filter specifications. This occurred 
when the program coefficients were reduced to three signif- 
icant octal digits in each case. Thus, eight magnitude bits 
or third-decimal-place accuracy was required for proper 
operation of the bandpass digital filter. 

The results of this example do not permit any general 
conclusions to be made regarding the choice of one filter 
form over another. It is likely that the arguments against 
the direct form from the point of view of coefficient accuracy 
which are made in the literature [1, 4, and 5] have greater 
significance in higher -order filters. Also the problem of 
filter stability was not encountered in this example, since 
the poles of the digital transfer function were located well 
within the unit circle in the z-plahe. 

Modification of the computer word length is accom- 
plished in the simulation program by inserting the ETR 
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instruction after each multiplication and altering the 
contents of location MASK as desired. Although the A/D 
converter provides only 13 bits of accuracy in the input 
data, all calculations are performed with 24-bit accuracy 
unless altered by the ETR instruction. Accordingly the 
effects of round-off error in a digital filter can be ob- 
served as the computer word length is reduced. 

The contents of location MASK were decremented from 
15 bits for each filter form. The point at which quantiza- 
tion noise reduces the response of the digital filter to an 
unsatisfactory state must be judged in light of the filter 
application. For this example eight bits and seven bits 
were required for the direct form and the parallel form 
respectively to retain basic signal definition on an oscil- 
loscope. The response of the cascade filter was unsatisfac 
tory even when the computer word was truncated to 15 bits. 
Figure 6-4 shows the output waveform for a 300 Hz signal 
using 15-bit arithmetic in the cascade and parallel filters 
Figure 6-5 shows the relative quantization noise in the 
spectra of the two signals. 

It appears that the cascade filter structure, with 
its serial form of signal processing, requires considerably 
greater accuracy in arithmetic computations than either of 
the other two structures. Hence, the significant advantage 
of using time-multiplexed circuitry is somewhat offset by 
the increased word length requirement. 
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300 Hz 





Figure 6-4. Quantization Effects in Digital Filters 
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C. A NOISE FILTER 



As a final example consider the problem of generating a 
noise signal whose power spectrum is of uniform height within 
a specified frequency range, but is significantly attenuated 
outside this frequency range. Such a noise signal has 
application in various signal-analysis and measurement schemes 
for communications systems, radar systems, etc. 

If a continuous-wave carrier signal is frequency-modulated 
by a gaussian noise source, the resultant signal spectrum is 
gaussian in shape and centered on the carrier frequency. The 
desired noise spectrum could be obtained if this gaussian- 
shaped spectrum could be filtered by a bandpass filter whose 
magnitude-ver sus-f requency characteristic varied as the 
reciprocal of a gaussian curve within the passband. 

Such a problem in spectral shaping can be solved using a 
digital filter. The bandpass digital filter which was pres- 
ented in Section III, C, as an example of matched z-transform 
design is offered as an approximate solution to the problem. 

The gaussian spectrum was generated by applying a low- 
frequency noise source {the amplitude distribution of the 
noise voltage is gaussian) to a voltage-controlled oscil- 
lator. The carrier frequency was chosen to be 450 Hz and 
the rms voltage of the noise was 1.0 volts. Applying a con- 
stant 1.0 volts to the VCO displaced the carrier frequency 
200 Hz. Therefore the gaussian spectrum centered at 450 Hz 
has a standard deviation of 200 Hz. 
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The bandpass second-order Chebyshev filter was 
designed for a 12-db gain between the peaks and the valley. 
This corresponded to a displacement of 1.6 standard deviations 
from the mean on the normal curve; or, in terms of frequency, 
for the given hardware it corresponds to the peak frequencies 
being separated from the carrier by 320 Hz. Thus, the ideal 
frequency characteristic for the specified conditions is 
shown in Fig. 6-6. The actual solution obtained by standard 




Figure 6-6. An Ideal Noise Filter. 

analog-filter design methods and a subsequent matched z- 
trans formation was plotted in Fig. 3-7. Notice that the s- 
domain solution in Fig. 3-6, which is itself distorted by the 
low-pass-to-bandpass transformation, suffers further distortion 
by the matched z-transformation. The final result is, however, 
a fair approximation of the desired frequency characteristic. 

The gaussian spectrum, centered at 450 Hz, is shown in 
Fig. 6-7a and the filtered spectrxam is shown in Fig. 6-7b. 
Notice that the gain factor in the digital filter was adjusted 
for 6db gain at the peak frequencies and a 6-db attenuation 
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a. Gaussian Input Spectrum 




b. Filtered Output Spectrum 



Figure 6-7. Output Signal and Spectrxam from a Digital 
Noise Filter, 
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at the center frequency. Further refinement in this filter 
design could be accomplished using optimization techniques 
in the s domain. 



VII . SUMMARY AND CONCLUSIONS 

The subject of digital signal processing continues to 
receive ever increasing attention. The basic appeal of 
digital circuitry is still lower cost, higher reliability 
and faster operation. Specifically, real-time digital filters 
now appear to be on the threshold of becoming a commercial 
reality. To a large extent the impetus for the theoretical 
development of the digital filter has been the promise of 
special-purpose computational hardware using LSI technology. 

It appears extremely likely that digital filters will find 
significant real-time applications in the not too distant 
future . 

An attempt has been made in this thesis to assemble and 
organize much of the current theory of digital filters, and 
to present in detail those aspects relating to the real-time 
implementation of digital filters. Furthermore, the proposed 
hybrid- computer simulation scheme provides a vehicle to 
process signals in real time and to observe the theory of 
digital filters in action. 

A. COMMENTS ON DESIGN METHOD AND IMPLEMENTATION TECHNIQUES 

Most of the theory of the z-transform and its application 
to discrete-time signals is not new. However, the development 
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of the digital filter has brought about new variations of 
that transform. For example, discussions of the bilinear 
and the matched z-transformi are not found in the older texts 
on sampled-data systems . Also the very popular digital- 
filter design technique of using s-plane continuous-filter 
synthesis has generated new and extremely useful "special 
purpose" transformations. One of these transformations, the 
bandpass version of the bilinear z-transform, obviates the 
conventional low-pass- to-bandpass frequency transformations 
in the s-plane. 

The usefulness of all the various transformations lies in 
their ability to be easily implemented and to produce accu- 
rate results. The use of large computer programs to synthe- 
size digital transfer functions from continuous transfer 
functions adds versatility to the indirect design methods. 

The three general implementation schemes which have been 
discussed in this paper all have useful applications. The 
direct form (including both the straight and canonical 
versions) is the simplest to use and gives very satisfactory 
results for relatively low-order filters. The cascade form 
represents a high-order filter as a cascade of low-order 
subfilters. This scheme is particularly adapted to realiza- 
tion in special-purpose hardware because of the "basic filter 
section" concept. A second-order filter section is time- 
multiplexed to permit the cascading of a single signal 
through the filter section n times. In this way a 2nth-order 
filter can be realized. 
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Time multiplexing is the key to the economical use of 
hardware. For example, the serial multiplier which is 
presented in Appendix A would seem to be a disproportionately 
complex circuit were it not for the built-in ability to 
handle multiplexed input signals and coefficients. 

The final filter structure is the parallel form. This 
form could provide the fastest signal processing if separate 
hardware were available for each subfilter. On the other 
hand, the parallel structure can also make use of time- 
multiplexed hardware. This implementation scheme retains 
the compact and economical construction of the cascade form, 
but could operate with a shorter word length. This is be- 
cause the cascade form inherently requires greater accuracy 
in the arithmetic unit. 

B. AN APPRAISAL OF THE REAL-TIME DIGITAL-FILTER SIMULATION 
The hybrid-computer program which was presented in this 
paper represents a first step toward the real-time processing 
of signals by digital filters. Theoretically any digital 
filter that is designed for real-time application could be 
simulated by this general scheme. (Real-time digital filters 
are considered to be of the recursive type in this context.) 

In the author's view there are two basic limitations to 
the digital-filter simulation scheme. The first is the 24- 
bit word length of the XDS 9300 computer- The coefficient 
and arithmetic accuracy required in high-order filters may 
not be attainable using a 24-bit computer word. A solution 
to this problem would be to rewrite the simulation programs 

for double-precision arithmetic. 
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The second limitation is that there is an upper limit 

■> 

on the allowable sampling frequency imposed by the cycle 
time of the XDS 9300 computer » Execution of a second-order 
difference equation (one basic filter section) in Computer 
Program 7 requires approximately 200 microseconds of computer 
time. This figure varies slightly depending on the arithme- 
tic sequence being used. This minimum execution time means 
that a fourth-order filter cannot be simulated using a 
sample frequency higher than 2500 Hz, or that an eighth- 
order filter cannot be simulated using a sample frequency 
higher than 1250 Hz, etc. From these figures it is evident 
that experimentation must be limited to the sub-audio range 
of frequencies. Indeed, the same general problem applies 
to all digital-filter systems. 

A partial solution to the problem at the Naval Post- 
graduate School is to write a more efficient program for 
the XDS 9300 computer. However, the future for the general 
simulation scheme and for digital filters as a whole is 
much brighter. The next generation of digital computers will 
make extensive use of LSI technology. Computers using mon- 
olithic circuits will have addition times in the order of 
50 nanoseconds vice 1.75 microseconds for the XDS 9300 com- 
puter. This reduction of almost two orders of magnitude in 
computing speed will mean an expansion of the digital- 
filtering spectrum to include the complete audio range and 
beyond. There is every reason to believe that special- 
purpose hardware will benefit from similar advances. 
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APPENDIX A 



A BINARY SERIAL MULTIPLIER FOR APPLICATIONS IN DIGITAL FILTERS 

1. INTRODUCTION 

Real-time implementation of digital filters involves the 
construction of a special-purpose digital computer. The 
arithmetic unit of this computer, which performs the multi- 
plication and addition of digital signals, presents some 
interesting challenges in logical design. One part of the 
arithmetic unit, the multiplier, is the subject of this 
Appendix . 

Specifically the computational aspects of an arithmetic 
unit which are peculiar to digital filters are discussed. 

From these considerations some general specifications for a 
multiplier are determined. Finally one solution to the prob- 
lem is presented using a binary serial multiplier. 

2. COMPUTATIONAL CONSIDERATIONS IN DIGITAL -FILTERS 

Most digital-filter transfer functions permit simultane- 
ous arithmetic operations during a given Nyquist interval. 

If the number of operations to be performed is large, economy 
can be achieved using serial arithmetic and by sharing the 
adders and multipliers. Time sharing in the filter can be 
accomplished with multiplexing circuitry. In addition to a 
significant simplicatibn of the hardware, serial arithmetic 
provides increased modularity and flexibility in the digital 
circuit configurations. Also the processing rate is limited 
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only by circuit speed and not by carry-propagation times in 
the adders and multipliers. Finally with serial arithmetic 
sample delays are realized simply as single-input single- 
output shift registers. 

The two's complement representation of binary numbers is 
also appropriate because addition can proceed starting with 
the least significant bit with no advance knowledge of the 
signs and with no later correction to sums being necessary. 
For example, the algebraic summation of the two signed num- 
bers, 1.101 and 0.010, can be accomplished as follows: 

conversion to 
two ' s complement 

1.101 0.011 

Subtract 0.010 ^ 1.010 Add 

0.011 1.101 

Addition using the two's complement form gives the correct 

signed result with no concern about signs or carries. In 

the multiplier, however, it is simpler to treat the sign 

separately and multiply only positive numbers. After the 

positive product is obtained the correct sign is reinserted 

in the proper bit location. 

Multiplication in digital filters is subject to a very 
severe restriction. Normally the multiplication of an N-bit 
number by a K-bit number results in a product (N+K) bits 
long. This situation cannot be tolerated in a serial system 
because the N-bit words (representing successive samples) are 
adjacent to each other in time. Consequently there is no 
room for the generation of a product in excess of N bits. 

If an N-bit sample is multiplied by a K-bit coefficient, the 
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product must be truncated and rounded to N bits. The solu- 
tion to this problem is to multiply one bit at a time using 
partial sums . It can be shown that the truncation of each 
partial sum and rounding the final partial sum will produce 
the same result as truncating and rounding a full final 
product . 

EXAMPLE 



Let X = the sample value = 0.00101, where the last bit 
in time is the sign bit (0 always) . 
a = the multiplying coefficient = 1.0110 with no 
sign bit. 

In this example N=6 and K=5. Then multiplication in the 
normal fashion can be compared to partial sum multiplication 
as follows: 
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The results are identically correct. 

If we require the magnitude of the sample to be less than 
one (i.e., |x|<1.0), and the magnitude of the multiplying 

coefficient to be less than two (i.e., |a|<2.0), then the 

fixed-point multiplication of the positive coefficient times 
the positive sample will always be of the form shown in the 
example. That is, the product is always positive and the 



111 



decimal point follows the sign bit just as ic did in the 
input to the multiplier. It is important to note that the 
sign bit of the coefficient is not included in the multi- 
plication. The restrictions on magnitude size are not 
unrealistic because the sample can easily be scaled. Also 
the multiplying coefficients for typical cascade and parallel 
filter forms are in the range -2.0 to 2.0. Any overflow in 
the partial sums will cause an error. This problem can be 
eliminated by further scaling and will not be discussed. 

In view of the above discussion the general requirements 
for the multiplier are the following: 

1. Accept an uncomplemented sample value from a serial 
delay device. This sample will be in binary form, N bits 
long, with the last bit in time representing the sign. The 
remainder of the word will be a fractional value less than 
one . 

2. Accept each single bit of a K-bit multiplying coef- 
ficient separately for use as a constant bit multiplier. The 
sign bit is used only to determine the sign of the product 
and the magnitude of the coefficient is less than 2.0. 



3. Remove the sign bit, from the sample value and 

multiply it by the coefficient sign bit, modulo 2 

(exclusive-OR) to determine the sign of the product. 



4 . Multiply the remaining positive sample by the positive 
multiplying coefficient using a serial multiplication scheme. 
The partial products must be truncated to N bits during each 
step and the, final N-bit product is rounded up by adding in 
the last truncated bit. The sign bit of this product must 
always be a zero (positive) . 



5. Reinsert the proper sign bit into the final product. The 
decimal point will be positioned immediately after the sign bit 
in the product. 
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6. Convert the signed product to the two's complement form 
for insertion into the following adder. 



3. THE MULTIPLIER 



A binary serial multiplier which will satisfy the above 
requirements is shown in block diagram form in Fig. A-1. An 
explanation of the multiplier follows: (Refer to correspond- 

ing step numbers on the diagram) . 

1. The notation of Fig. A-1 is based on a sample value N 
bits long including sign and a multiplying coefficient K bits 
long not including sign. Time "t^" is defined as the time the 
first bit of the sample, , becomes available at the input. 
This is the least significant bit of the sample word. Then 
succeeding times become t , ^2 , * * * tj^_ (which is the time the 
last bit of the sample, x^^, becomes available at the same 
point) . 



2. First the sign bxt of the sample is removed and that 
bit is converted to zero (positive) in the sample word. This 
operation normally requires an N-bit delay but the delay will 
be disregarded for sake of clarity in the discussion which 
follows . 



3. The removed sign bit, x^^, is multiplied modulo 2 with 
the sign bit of the coefficient, • The resulting product 

sign is stored or delayed for later use. 



4. At time t^ the least significant bit of the sample, x^ , 
is multiplied by the least significant bit of the coefficient, 
aj^. This product is the first bit of the first partial sum. 
The inputs to the multiplier bit sections labeled (j=0,l, 
•••K-2) are timing signals which accomplish the truncation. 

At time t^, r^ would be low (the other being high) so the 
output of partial sum 1 at time t^ would be gated to zero. 

At time tj^ r^^ goes low so the output of partial sum 2 at t^^ 
is zero. In this fashion the signal at the input to the 
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Multiplier Bit Sections 





Note: 



r j is low at bit time j . The line is high all other 
times . 



Figure A-1. A Binary Serial Multiplier. 
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multiplier follows the sequence Xj^(t^), X 2 ( t^) , • • *Xj^ ( 
while the signal at the output follows the sequence O(t^), 

0 (t^) , • • *0 (tj^_2> / ^2 * *^N ^^K+N-2^ ' 

illustrate the scheme more clearly consider the following 

example : 

EXAMPLE 

Let x=0.010 N=4 

a=1.011 K=4 

Then x = 0.010 
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(K-1) bit delay in product. 

It is easily seen that the preceding sample, occurring at 
the input times , t_^, ^-1 have produced a 

correct product at the output at times t^, t^^, t 2 . 

5. In the last multiplier bit section, the bit which is 

truncated at bit time t^ „ is fed back through the carry 

J\“” 

loop to accomplish the rounded product. 

6. Finally the sign bit (zero in all cases except over- 
flow) is changed to the correct value and the correctly 
signed product is then converted to two's complement form for 
insertion into an adder - 

It is seen from Fig. A-1 and the above discussion that 
this serial multiplier has an inherent delay of (K-1) bits, 
not counting any delay encountered in the sign bit removal. 
This additional delay must be deducted from the prescribed 
delay in the multiplication branch for proper operation of 
the digital filter. 
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OOOOOOOOOOO 



COMPUTER PROGRAM 1 



EREOUCMCY RESPONSE OF A DIGITAL TRANSFER FUNCTION 

M IS THE NUMBER OF CASCADED SECOND-ORDER TRANSFER 
FUNC TIONS 

FI AND F2 ARE THE LOW AND HIGH FREQUENCIES IN THE 
FREQUENCY RANGE OF INTEREST 
T IS THE SAMPLING INTERVAL 

MAGH IS THE MAGNITUDE OF THE TRANSFER FUNCTION 
EVALUATED AT Z = EXP ( - J*OMEGA=»=T ) 

MAGNH IS MAGH EXPRESSED IN DECIBELS 

DIMENSION A00(10),A11(10),A22(1C),B11(10),B22(10) 

REAL MAGH, MAGNH, NUMl ,NUM2 
INTEGER F1,F2 
COMPLEX NUM,DENOM 
RFAD(5,?CC)M,F1,F2,T 
200 FORMAT! 12, ?IA,F10.0) 

WRITE(6 ,1CC) 

ICO FORMAT ( IHl , lOX , 'FREQUENCY' , 5X, • ABS GA I N ' , 6 X , ' G AI N ( DB ) ' 
^/) 

DC 5 1=1 ,M, 

READ(5, 101 )A00 ( I) , All ( I ) , A22( n ,B1 1 ( I ) ,B22 ( I ) 

101 FORMAT (5F1C,0) 

5 CONTINUE 
F = F1 

DO 1C I=F1,F2 
X=F^6. 2821 85306-*T 
MAGH=1.0 
DC 6 J=1,M 
AC=AOO( J) 

A1 = A11 ( J ) 

A2=A22( J) 

B1=B11( J) 

B?=B22(J) 

NUM1 = ( A2 + AO*COS(X )+Al 
NUM2=(A0-A2)’*‘SIN(X) 

NUMo=CMPLX( NUMl ,NUM2) 

DEN0M1=(P2 + 1« )*'CDS(X)+B1 
DEN0M2=( l.-B2)*SIN(X) 

DFNnM=CMPlX(DENOMl ,DEN0M2 ) 
MAGH=APS(MAGH*CABS(MUM/DENOM) ) 

6 CONTINUE 

MAGNH=2C,» ALOGIO(MAGH) 

WPITE(6,1C2) F,MAGH,MAGNH 

102 FCRMATdH , 09X , F6, 1 , 8X , F8 . 6 , 7X, F g. 3 ) 

F=F+1.C 

10 CONTINUE 
STOP 
END 
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COMPUTER PROGRAM 2 
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COMPUTER PROGRAM 2 (Continued) 



UJ 


















« — 1 












+ 


:z 


















1— UJ 






Ul 




UJ 


r>Qc: 




UI 


>r 




_j 


CLC 




3: 






CL 


»-CL 


+ 


»— 1 






31 


Z>UJ 


z 








< 


<DCr 






cr 




O) 


CL 


UJ 


Ui 


(D 






u. 


z: 


_j 


u. 




or 


(D(D 


»— • 


CL 






<D 






3" 


in 




u. 


< 




< 


CL 


o 




\ < 


Ui 


in 


OD 


\ 


3T 


Q»— 


-j 






< 




C 


GL 


(Y 


lU 




CO 




z: 


<D 


u 


z: 




cr 


< 


U» 


z 


cr 


_J 


<2DI— 


cn 




UJ 


CD 


c 


b-U 




>- 


X 


U. 


2: 




cr 




r: 


cr 


» — « 


UiX 


(2D 


< 


(D 


UJ 


u. 


CL CO 


u. 




U 


CL 



OJ 

< 1 - 

o 

O OJ 

o 











o 


CVJ ^ 




CD CD 




o 










a« 




(or^ 




CVJ 


-J-J 


-- CD 


o z 


CD 


CO 




-J* 






CLrv 


< < 


LO 


U.Ll 


LO _l 


O CD 


_J 


CL 


T 


3“ Z 


z 


z 


cnr^ 


1- »- 


CVl 




u 


^ C J 


U- 


cn 


(Vj -D 




"DOJ 




<o 


< < 


1 o 


_J_I • CD 


o o 


on o 


O 1 


o 


< ^ (/) 


CO CO ^ CO 


CO cn 


•*— t CO 


O M 


o n 




XX CA (C 


— c 


O C 


< 40 


<c 



n 

_i (O z 


c J 


o 

CO 3’ 


< J 


O 

COZ 


Z X 


C < 


z 
X CL 


c z 


IDO 


UJ 


a <f 


Z K 


Z 3 


X 


ID JO 


oo 


_IQ 


O ID 


_JQ 


xo 


OJ- 


X CD 


»— X 


XX 


fNl 


CD 1— 


CD CD 


X 


X 


> < <t 


_J> 


< < 


_J3 


< < 


CQJ 


_JC0 


CD O 


COCO 


CD m 


a. 


oco 


Lua 


CO CD 





ojon 


^ LO 


vo 


00 (Tv 


O-H 


cvion 


^ in 


nO r- 


00 (T) 


o-.-^ 


CO 

a 

ID 

e 

< 

CVJ on 


4- in 


vor- 


00 o> 


♦ 

O 


on onm 


on on 


on on 


on on 




J 




J 




mio 


UDLO 


in in 


in in 


Lfj in 


vO nC 





1 


vD 




CD 




nO 


O 


1 — < 


vO 




4 




X 




o 


in 


in 


4 


sO 


o 


o 


4 


CD 


cn 


4 


4 




on 


o 


X 


X 


4- 


o 


OJ 


on 


4 


CDX 


vO 


4 


m 


m 


4 




X 


X 


in 


^ ^ 


o 




X 


O 




X 


vO 


in 




o 






^ — k 


CD 


T-W 


« — 1 


^ — i 


o ^ 


o 






^ — i 


o 








o 


o 


o 


r- 


f 


c 




tH 


o 


o 


o 




o 


o 


O 




o 


o 


o 


o 


o 


O 


o 


o 


o 




o 


o 


o 


o 


o 


r^ 


O 


4 


o 


C 


o 


o 


o 


4" 


CD 


o 


o 


4 


o 


o 


o 


4 O 


o 


CD 


o 


o 


o 


X 


o 


CD 


o 


CD 


o 


cn 


O 


X 


CJ 


o 


o 


o 




o 








O 








O 






































on 


o 


IT 


nO 


on 


O 


IT 


vC 


on 


oin 


m 




vO 


nC 




o 


vC 


X 


•w-H 


^ — i 


o 


o 


vC 




n — f 


X 


f-l 


T~i 


vO 


vO 


on 


i 


vO 


nO 


on 




vO 


vom 


o 


w-i 


^ — 1 


rv 


in 






in 


o 


O 


o 








X 


in 


o 


4 
































o 












o 




CVj 










O 


o 


o 


o 


o 


o 


o 


o 


o 


oo 


ox 


X 


X 


X 


4 


o 


o 


o 


O 


o 


4 


o 


o 


CD 


o 


CD 


CD 


O 


CD 


o 


o 


CD 


o 


o 


o 


CD 


oo 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


o 


O 


o 


o 


o 


o 


o 


o 


a 


on 


4“ 


in 


vO 




c^ 


^ — i 


X 


004 


in 


vO 




o 


1 


X 


X 


4 


UD 


vO 




CD 


1 — < 


X 


X 


4 


U') 


vO 


on 


on 


X 


m 


on 


on 


4 


4 


4 


4 4 


4 


4 


4 


in 


in 


in 


LO 


in 


in 


in 


in 


vC» 


xO 


vO 


vO 


vO 


vO 


vC^ 


o 


o 


o 


o 


CD 


CD 


O 


o 


O 


OO 


O 


o 


O 


o 


o 


o 


o 


o 


o 


o 


o 


CD 


CD 


C ) 


CD 


o 


CD 


o 


o 


o 


o 


o 


O 


o 


O 


CD 


o 


OO 


o 


O 


o 


o 


o 


o 


O 


o 


o 


o 


o 


O 


o 


o 


o 


o 


O 


CD 


o 


o 


o 


o 


O 


o 


O 


o 


o 


oo 


O 


O 


o 


o 


o 


o 


o 


o 


o 


o 


o 


O 


o 


CD 


CD 


CD 


o 


o 



118 



Q0067 0 0 00 00000 62 DA0PS PZE PERF0RM D/A 

00C70 0 C 16 C0126 63 LDA SUM 



COMPUTER PROGRAM 2 (Continued) 



o 

<D 



UJ 








h- 




CG QC 




UJ 




UJ 








or CL 








<: 






< 


UJ 




<n 




CG 


(n 


u 


0 


0. 








_j x: 


CL 


<D 


a: 


< < 


3 


K- 


(D 


z tn 


tr 




Lu 


CD 


cr 


O 




— j- 


UJ 


\ 


h- LU 


{/} X 


h- 


< 


•— 1 — 


Ul 






< LJ 


v: 2 T 




cr 


3: _J 


CJ 




<D 


G- 


CD cr 


Ul 




5: 


^ CD 


_J 




<D <E> 


CJ u. 


CD 


♦- 


7 CJ 




<C 


t— • 




: 3 c CG 


CO 


< 


(D CD 


Ul CL 


►— i 


3: 


0 h- 


^ CD 


0 











0 




X 


X 






0 


0 


Q 


CD 


z 




CD C 


CD 




0 


tr 5: 0 z CO 


CD 


2 


_J 


UJ 




CJ h- 


u 


C 3 Q 


X 


V 0 <D 0 (D a 


_J UJ 


Ul 


u 


h- 




in 0 < 


c 




Ul 


(o <: CJ tn CJ CD 


Ll ^ n 


c 


1 — ^ 


-J 


1 


'•-» < 0 


n 


UlLil 




< c <r fn -< < 


0 1 < 


0 


_J 1 


X 


<n 


^ ^ 


•K 


o< 




X 0 0 0 0 0 


* 


* 


X 


★ 


(D 


C^ ^ 


r-i r — i 


<ro 


T 



j :. 



or 


0 < 


X 1 - 


cr 


Ul cr 


3 CJ 




LiJ u 


Ui 


cr 3 


CJ 


cr 3 




or 


I?' 


UIZ 


(G UJ 


5 >■ 


X 


►- 


0 •- 


cD cn 


cr 


M V 


DC 


cr 




M cr 


IS! 


V cr 


cr 


— *cr 




CD 


CD CD 


M(D 


UJM 


cr cr 


cr 


u 


< (O 


U) CL 


rc 


Cl (O 


CD 


CD 




CL CD 


a 


C 0 X‘ 


cc 


n CD 




U 


CJ CJ 


CL CJ 


cr CL 


X'CT'' 


rn 




















0 






7 














H- 










n 








0 








3 






Z X 


z 


X 


Q< 


3 










y 










LjU 






or 






Cv (t; 


CD 


CD 


<0 


T 










UJ 








UJ 


h- 






0 






CJ CJ 


CJ 


U 


XX 


>' 










C J 








< 


_l 






X' 




CT; 


0 0 


< 


< 


cr cr 


cr 








♦ 


< 






Hf- 


0 


♦ X 




♦ 


UJ 


* 


CJ 


< c 


0 


Oi 


mm 


m 




ID vO 


1 ^ 00 


Ch 0 


T-l OJ 


cn 


^ ID 


vO 


00 U\ 


0 «-l 


c\jcn 


^ in 


vO 


oc 0^ 


O-rH 


ojcn 


^in 


vO 


xO 


0 sO 


vC vO 


vo r- 


rv 












00 00 


00 00 


oc c)0 


oc 


00 


OG 00 


cno> 


cn (n 


Chcn. 


cn 



Co 




0 ^ 0 * 


Ov 


CJ ^ 


0 


n- 


0 m 


0 


UJ vO 


UJ 


CO 




0 


0 


PvP) 


in 




CG ^ 


0 


xO 


CJ a 0 


r- 


0 0 


0 


(VI 0 


0 






0 


0 


r^o 


0 






CJ 


CJ 






0 


0 ^ 


0 


»— < ■» — 1 




0 




0 


0 


0^ 


^ 1 — ♦ 


0 


0 


IP 0 


0 


0 0 


0 


0 


0 CJ 


0 


00 


0 


0 


4- ^ 


0 


0 


CJO 


0 


0 


0 u 


c^> 0 


0 


CJ 0 


0 


0 


C) 0 


0 


OCJ 


CJ 


^ 0 
0 


^tH 


C J ^ 
^—1 


c: 


CJO 


CJ 




Li' xC 


O'* ^ 




0 m 




r- 


0 


0 


rn ^ 


1^ 


0 -H 


r.' c* 


c-c 


0 


pp 


p 




0 


m 

fVJ 




0 


0 


UJ 


0 UJ 


0 


P-O 


10 


00 

a 


0 0 


00 


0 


OCJ 


0 


0 


CJ CJ 


0 0 


0 


0 

0 


CJ 


CJ 


0 0 


0 


0 

0 


0 


CUCJ 

0 


' -4' 

0 0 


CJ 

0 


0 


00 


C.J 


CJ 


0 0 


0 

0 


0 


0 

0 

0 


▼H 


0 ^ 


0 


00 


«— * 


00 


0 0 


00 


0 


00 


0 


-r— • 


a; a 


3 UJ 


vO 


0 




a 




UJ 


xo 


CJ 


▼-» a? 


cn 4- 


UJxXJ 


rv. 0 




p 




rv 


rv. 




0 


0 


0 


CJ 0 


0 


00 




vH ^ — < 


^ — ! 


^•1 1 




a a* 


a 


0 


CJ 0 


0 0 


0 


0 i 




— 1 


^ 1 


▼H 




T-» 


V— # 


^ # 


■* — 4 








CJ 


CJ CJ 


CJ C' 


O' 


0 CJ 


0 


CJ 


0 0 


0 


0 0 


0 


0 (^ 


0 0 


00 


0 0 


OCJ 


CJ 


0 


CJ 0 


c 0 


CJ 


0 0 


0 


0 


CJ CJ 


CJ> 


0 0 


0 


0 CJ 


0 0 


00 


0 CJ 


OCJ 


0 



119 



COMPUTER PROGRAM 2 (Continued) 



c 



x: 


1— 


z> 


< 


in 


Q 


(X 


or 


(b 


<D 


u_ 


li_ 


Ul 


bJ 


<3 


O 


< 


< 


q: 


(T 


(D 


O 


1— 


l-CO 


o> 


o)(n 




UJ 


UJ 


uitr 


1— 


1—0 


c 


C Q 


u 


u<c 


<D 


(t) 


_J 


_JC 


_J 


_J\ 


c 


<c o 







o 








mo 






vO 


^o 






r^o 


4-rv 




cvjin 


OJO 




h- 


(n 00 


Ch ^ 


rorv 


cr 


cum 


cuo^ 




<■ 


00 >o 


00 


-»-• 


•— 


1 

1 

3 

3 

5 

2 

5 


(U ^ 


» o 


cn 







c 


c < 


<r c 


cc 




01 


men 


cni— 


1— >- 


►~h“ 


h- 1— 


o 


UJ 


UJUl 


Ul< 


c <c 


c c 


c c 


z 


or 


ct ct 


(X o 


oo 


oo 


oo 


UJ 




CD 












o 


-J 


ex 










-j 


u_ 


<t Q 






V 




li- 


•- x: 


1— 






CO 




o 


_iZ) 


c c 


o ^ 


(U^ 


(U c 




c 


xen 


oo 


< < 


CCD 


mx 






00 cn 


o^ 


OJ cn 


-:rm 


Nor^ 


00 


CTN 


cnen 


oo 


oo 


oo 


oo 


o 








«-» ^ 


^t4 ■•—I 


f 1 


w-t 



lo vo no 
o in ro LDOj 
o ^ ^ OJ 
o o ^ o m 
o m (V! in oj 

O O ^ O (VI 
O ^ <\i ^ CVj 
O O O O OJ 



D o o in ^ 
mo o o o 
OJ o o o o 
ojr^ o o o 
-r-lfV o o o 
o rv 
o ^- 
^v o r^ 



-^mvo«r-i^mvor^o^(Vj 
ojcvjojmnroncn^^^' ^ ^ 



ooooooooooo oo 

ooooooooooo oo 



120 



COMPUTER PROGRAM 3 



1 + DIGITAL FILTER SIMULATION 



2 ♦ SFCONO-ORDER LOW-PASS 


3 ♦ 






♦ 












♦ 






* 












21 


BRM 


ADOPS 


22 


LDA 


SUM 


23 


APSA 


2 


24 


STA 


SUM 


25 


LDA 


SUM + 1 


26 


MUL 


B1 


27 


ALSO 
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28 


ADM 


SUM 


29 


LDA 


SUM+2 


30 


MUL 


B2 


31 


ALSO 
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32 


ADM 


SUM 


33 


LDA 


SUM 


34 


MUL 


AO 


35 


ALSO 
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36 


STA 


Y 


37 


LDA 


SUM+1 


78 


MUL 


A1 


31 


At sn 


1 


40 


ADM 


Y 


41 


LDA 


SUM+2 


42 


MUL 


A2 


43 


ALSO 


1 


44 


ADM 


Y 


45 


BRM 


DAOPS 


46 


LDX 


=077700002 


47 


LDA 


SUM, 2 


43 


STA 


SUM+1 ,2 


49 


BRX 


$-2,2 


♦ 












ilc 

-V 






-JC 






*e 






* 






w- 






♦ 






102 SUM 


RES 


3 


103 Y 


PES 


1 


104 DAACR 


DATA 


5 


105 AO 


DATA 


282927 


106 A1 


DATA 


565854 


107 A2 


DATA 


282927 


1C3 81 


DATA 


4794006 


109 B2 


DATA 


-1731413 


no MASK 


DATA 


077777000 


111 


END 


START 



FILTER (CANONICAL FORM! 



TAKE SAMPLE 
DIVIDE SAMPLE BY A 



OUTPUT VALUE AT TIME N 
PERFORM 0/A OF OUTPUT 
SHIFT DATA TO PREPARE 
FOR SAMPLE TIME M+1 



0/A ADDRESS 
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31 
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39 

40 
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45 

46 

47 

48 

49 

50 

51 

52 

53 
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53 
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60 
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117 

113 

119 
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121 
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127 
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COMPUTER PROGRAM 4 



CIGITAL FILTER SIMULATION 
* FOUPTH-ORDFR BANDPASS FILTER (DIRECT FORM! 








PPM 


ADOPS 




LDA 


DATA 




ARSA 


2 




sta 


DATA 




MUL 


AO 




ALDS 


1 




STA 


SUM 




LDA 


DATA+1 




MUL 


A1 




ALSO 


1 




ADM 


SUM 




LDA 


DATA+2 




MUL 


A2 




ALSO 


1 




ADM 


SUM 




LDA 


DATA+3 




MUL 


A3 




ALSO 


1 




ADM 


SUM 




LDA 


DATA+4 




MUL 


A4 




ALSO 


1 




ADM 


SUM 




LOA 


SUM+l 




MUL 


B1 




ALSO 


1 




ADM 


SUM 




LDA 


SUM + 2 




MUL 


B2 




ALSO 


1 




ADM 


SUM 




LDA 


SUM+3 




MUL 


B3 




ALSO 


1 




ADM 


SUM 




1 DA 


SUM + 4 




MUL 


B4 




ALSO 


1 




ADM 


SUM 




BPM 


DAOPS 




LDX 


=077700010,2 




LDA 


data, 2 




STA 


DATA+1 ,2 




ERX 


$-2,2 


» 






lie 












CAT A 


RES 


5 


SUM 


RES 


5 


DAACR 


DATA 


5 


AO 


DATA 


549898 


A1 


DATA 


000000000 


A2 


DATA 


-1099797 


A3 


DATA 


000000000 


A4 


DATA 


549898 


B1 


DATA 


5872294 


B2 


DATA 


-5336061 


B3 


DATA 


2761609 


B4 


DATA 


-1141761 


MASK 


DATA 


077777000 




END 


START 



TAKE SAMPLE 

DIVIDE SAMPLE BY 4 
DATA IS NOW THE INPUT 



SUM IS ACCUMULATED OUTPUT 
AS TERNS OF DIFFERENCE 
FQN ARE ADDED TOGETHER 



OUTPUT VALUE FOR TIME N 
PERFORM D/A OF OUTPUT 
SHIFT DATA TO PREPARE 
FOR SAMPLE TIME N+1 



D/A ADDRESS 
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COMPUTER PROGRAM 5 



1 ■* DIGITAL FILTER SIMULATION 

2 * FOURTH-OPOER BANDPASS FILTER (CASCADE FORM) 



3 






* 






* 












♦ 






* 






21 


BPM 


ADOPS 


22 


LDA 


SUM 


23 


ARSA 


2 


24 


STA 


SUM+1 


25 


LDX 


=077700001 


26 


COPY 


{0»5) 


27 REPEAT 


STA 


SUM 


28 


LDA 


SUM+2, 3 


29 


MUL 


B1 f3 


30 


ALSO 


1 


31 


ADM 


SUM, 3 


32 


LDA 


SUM+4,3 


33 


MUL 


B2,3 


34 


ALSD 


1 


35 


ADM 


SUM, 3 


36 


LDA 


SUM, 3 


37 


MUL 


AO, 3 


38 


ALSD 


1 


39 


STA 


Y,3 


40 


LDA 


SUM+2, 3 


41 


MUL 


A1 ,3 


42 


ALSO 


1 


43 


ADM 


Y,3 


44 


LDA 


SUM+4,3 


45 


MUL 


A2,3 


46 


ALSD 


1 


47 


ADM 


Y,3 


48 


LDA 


Y,3 


49 


BPX 


REPEAT, 3 


50 


BPM 


DAOPS 


51 


LDX 


=077600004 


52 


LDA 


SUM, 2 


53 


STA 


SUM+2, 2 


54 


LDA 


SUM+1,2 


56 


BPX 


S-4,2 


55 


STA 


SUM+3, 2 


* 






* 






♦ 


















A 






A 






* 






109 SUM 


RES 


6 


110 Y 


RES 


2 


111 DAADP 


DATA 


5 


112 AO 


DATA 


1518695 


113 


DATA 


1518695 


114 A1 


DATA 


3037389 


115 


DATA 


-3037389 


116 A2 


DATA 


1518695 


117 


DATA 


1518695 


118 B1 


DATA 


746758 


119 


DATA 


5125536 


120 B2 


DATA 


-1890880 


121 


DATA 


-2532626 


122 MASK 


DATA 


077777000 


123 


END 


START 



/ 

TAKE SAMPLE 

DIVIDE INPUT VALUE BY 4 



PERFORM D/A OF OUTPUT 
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COMPUTER PROGRAM 6 



EVALUATION OF RESIDUES IN A PARTIAL FRACTION EXPANSION 

N IS THE ORDER OF THE DENOMINATOR POLYNOMIAL 

A(I) ARE THE REAL NUMERATOR COEFFICIENTS IN ASCENDING 

ORDER BEGINING WITH THE CONSTANT TERM 

Pin APE THE COMPLEX POLE FACTORS. THE MULTIPLICITY OF 
POLES CANNOT EXCEED 2 

K IS THE COMPLEX RESIDUE ASSOCIATED WITH A GIVEN POLE 
FACTOR AS COMPUTED BY SUBROUTINE RES 

COMPLEX K,P 

DIMENSION A(6) ,P( 6) ,K( 6) 

1 READ(5,100)N 

100 FCRMAT(Il) 

IFIN.EQ.OGO TO 99 
READ(5tlGl )(A( I)tl = l,6) 

101 FCRMAT(6F10.0) 

READ(5, 1C2) (P( I ) , 1=1,6) 

102 FCRMAT(8F1C.0/AF10.0) 

CALL RFS(N,A,P,K) 

WRITF(6,20C)N, (A( I ) ,1=1,6) 

200 FCPMATdHl ,6X,*THE DEGREE OF THE DENOMINATOR 
’^POLYNIMIAL IS •,I1//6X,*THE COEFFICIENTS OF THE 
*NUMERATCP POLYNOMIAL AR E: • / / ( 13X, 6F 10. 6 ) ) 

WRITE(6,2C1) (P( I ) ,K(I ) , I=1,N) 

201 FORMAT! IHC, 16X, 'POLE F ACTOR * , 10 X , • R E S IDUE S* //1 5X , 

= *RrAL* ,6X, • IMAG*,6X,*REAL*,6X, • IMAG*//( 11X,4F10.6)) 

GO TO 1 
99 STOP 
END 



SUBROUTINE RES(N,A,P,K) 

COMPLEX K,P,S, Z,NUM,DENOM,FUNCS,FUNCDS, DELTA 
DIMENSION A(6) ,P(6) ,K(6) 

NUM( Z)=A1+A2*Z+A3*Z**2+A4*Z**3+A5*ZY*4+A6*Z**5 
A1=A(1) 

A2=A(2) 

A3=A(? ) 

A4=A<4) 

A5=A(5) 

A6=A(6 ) 

DELTA=( .001 ,0. ) 

1 = 1 

5 L=0 

6 S=-P(l ) 

7 DENOM=( l.,0. ) 
no 50 J=1,N 

IF(CAPS(S+P( J) )-.05)9,9,8 

8 DENOM=DENOM*(S+P( J ) ) 

GC TO 5C 

9 L=L+1 

50 CONTINUE 

GC T0( 11,12,13,13, 14, 14), L 

11 K(I)=NUM(S)/DENOM 
GO TO 10 

12 K( I ) =NUM(S) /DENOM 
1=1 + 1 

GC TO 6 

13 FUNCS=NUM( S) /DENOM 
S=S+DELTA 

GO TO 7 

14 FUNCDS=NUM(S)/DENOM 

K(I ) = (FUNCDS-FUNCS)/DELTA 
10 IF( I.EO.N) GO TO 90 
1 = 1 + 1 
GO TO 5 
90 RETURN 
END 
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COMPUTER PROGRAM 7 



1 * DIGITAL FILTER SIMULATION 

2 * FOURTH-ORDER BANDPASS FILTER (PARALLEL FORM) 

3 * 

* 

4c 

♦ 



21 


BRM 


22 


LDA 


23 


ARSA 


24 


STA 


25 


MUL 


26 


AL SD 


27 


STA 


28 


LDX 


29 REPEAT 


LDA 


30 


STA 


31 


LDA 


32 


MUL 


33 


ALSO 


34 


ADM 


35 


LDA 


36 


MUL 


37 


ALSO 


38 


ADM 


39 


LDA 


40 


MUL 


41 


ALSO 


42 


STA 


43 


LDA 


44 


MUL 


45 


ALSO 


46 


ADM 


47 


LDA 


48 


BRX 


49 


LDA 


50 


ADM 


51 


LDA 


52 


ADM 


53 


LDA 


54 


BRM 



ADOPS 

DATA 

2 

DATA 

GAMMA 

1 

ANSWER 

=077700001»3 

DATA 

SUM, 3 

SUM+2,3 

Bl,3 

I 

SUM, 3 
SUM+A, 3 
B2,3 

I 

SUM, 3 
SUM+2,3 
A1 ,3 
I 

Y,3 

SUM+A,3 

A2,3 

I 

Y,3 

Y,3 

REPEAT, 3 

Y 

ANSWER 

Y + 1 

ANSWER 

ANSWER 

ADOPS 



TAKE SAMPLE 
DIVIDE INPUT BY 4 



PERFORM 0/A OF OUTPUT 



4c 

♦ 

♦ 

4c 

4c 



X 

* 



4c 

4c 

4c 



13 


DATA 


PFS 


1 


14 


SUM 


RES 


6 


15 


Y 


RES 


2 


16 


ANSWER 


RES 


1 


17 


DAAOR 


DATA 


5 


18 


GAMMA 


DATA 


549898 


19 


A1 


DATA 


-1609581 


20 




DATA 


2379471 


21 


A2 


DATA 


627585 


22 




DATA 


194230 


23 


B1 


DATA 


5125536 


24 




DATA 


746758 


25 


B2 


DATA 


-2532626 


26 




DATA 


-1890880 


27 


MASK 


DATA 


077777000 


28 




END 


START 



D/A ADDRESS 
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